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This thesis discusses the modif icaticr. of an e x i st ing 
Cptimal Systems Conttcl FOR IRAN Program (OPTSYS) originally 
obtained frcm Professor Arthur E. Bryson of Stanford 
University - 

The notified FORTRAN program (0PT5YSX) is now designed 
to run completely interactively under Vil/CMS on the IBM 3C33 
and is considered completely compatible with similar oper- 
ating systems. 

Program capabilities include: complete eigensystem 
analysis; the ability to perform computations on very large 
multivariable systems; controller, filter, regulator and 
compensator synthesis; transfer function analysis; steady- 
state gain determination; and medal analysis. 

The program permits users to rapidly carry out simula- 
tion, analysis, and design cf all classes of Optimal Systems 
Centre! problems in a totally interactive mode. Examples of 
various types of problems are worked out during individual 
terminal sessions. 
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SYMBOLS 



A 

g 

C 

D 

E 

F 

G 

GAM 

E 

K 

he 

Kg 

Ns 

So 

F 

C 

E 

S 

u 

V 

w 

vO 

X 

x e 

y 



state (Ns, Ns) or output (No, No) weighting matrix 
control (lie, No) weighting matrix 
control gain matrix ( N c , N s ) 

control (N c ,Nc) or ncise (No,Ng) feedforward 
matri x 

expected value 

cpen-lcop dynamics matrix (Ns, Ns) 

control distribution matrix (Ns,Nc) 

state disturbance distribution matrix (Ns,Ng) 

measurement scaling matrix (No, Ns) 

estimator gain matrix (Ns, Nob) 

number of controls 

number of process ncise sources 

number of states 

number of observations or measurements 
covariance matrix of estimate error (Ns, Ns) 
white process noise covariance matrix (Ng,Ng) 
white meas. noise covariance matrix (No, No) 
steady-state covariance matrix of control (Nc,Nc) 
control vector (Nc,1) 

white measurement ncise vector (No, 1) , with 
zero mean and covariance matrix R 
white process noise vector (Ng,1), with 
zero mean ar.d covariance matrix Q 
constant disturbance vector (Ng, 1) 
state vector (Ns, 1) 
estimate of state vector (Ns,1) 
c utp ut/mea suremen t vector (No,1) 
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:n tbcduction 



1 be purpose cf -bis ohesis is to describe the extensive 
nodi fic a tier and improvement of an existing FORTRAN program 
(OPT SYS) used in the study, design, and application cf 
Optimal Systems Control theory. 

This optimal systems control program was originally 
developed by Hall [Eef. 1] in 1971 to support his research 
in rotary-wing VTOL aircraft control systems. The latest 
program modifications were made by Walker [Ref. 2] and Liu 
[Ref. 3] cf Stanford University, to OPTSYS 4 and CFTSYS 5 
respectively. These program versions performed quite satis- 
factorily in the batch environment, but exhibited varying 
degrees cf user hostility due to data input format require- 
ments and incomplete program documentation. 

. The original intent of this work was to adapt Walker’s 
modified version of CPTSYS to run interactively under VJS/CHS 
on the IEM 3033; however, the extensive modifications of 
CPTSYSX row represent a high-level interactive applications 
software system capable of integrated simulation, analysis, 
synthesis ar.d design of broad classes of optimal systems 
control problems. With OPTSYSX users may now evaluate 
various specialized optimal systems control applications, 
relieved of the burden of lengthy mathematical program 
development . 

It is assumed that the reader/user is familiar with the 
basic ccrcepts of Ccrtrol Theory and Optimal Systems Design. 
The problem descriptions and program development fellow the 
terminology and symbol/ naming conventions of Erysot 
[Ref. «]. 

An explanation cf the basic system of equations, the 
terms and symbology used, and a program overview including 
the general methods cf solution are presented first. 
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Interactive program dev sic c meat i 






data ::pu*. 

Inis *crk concludes with eranpie 
problems deicn stratec in the interac 
copy cf each terminal session with 
complete program listing is included 



s then discussed, with 
options available for 




tive mode, including a 
the final results. A 
in Appendix A. 
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II. THE QPTSYSX COMPOTES PROGRAM 
A- INTECEQCTION 

CP1SYSX is a double-pre cision FORTRAN program employing 
modern control theory analysis techniques. Although the 
program (.as originally written to syntnesrze controllers for 
rotary-wing VIOL aircraft [ 5ef . 5], it has been extensively 
modified tc enable controller, filter, and regulator 
synthesis as well as transfer function and modal analysis on 
ether types of large, mul ti- variable systems cf equations. 
The program modifications described in this work new allow 
rapid numerical computer analysis in a completely interac- 
tive mode. 

E. S3SIZM/K0DEL DESCRIPTION 

CEISYSX treats a linear stationary system model: 

x = [F]x + [G]u + £GAMMA]w (2.1) 

output equation 

y = £ H ]x + [D]u (2.2) 

measurement equation 

z = [HjX + £ D ] w + v (2.3) 

where 

u = ccntrcl vector (m X 1) 

v = white measurement noise vector (p X 1) 
t = white process noise vector (q X 1) 
state vectcr (n X 1) 



x = 



y = output vector (p X 1) 
z = measurement vector {p X ij 

r r * • rr 4. c* ^ ^ c ' »-« <-v v~ -1 v *■' t ~j ^ £ c: *•») .- — r~ *? w / c; *7 c ■*" "n r 7 — - - ^ ~s q — 

plan* matrix); [G] is the control distribution matrix; 
[GAMhA] is tha process noise distribution matrix; [H] is the 
measurement distribution matrix; and and [D] nay represent a 



feed-fc 


r war 


c distr 


ibution m 


atrix cf either 


the process noise 


vector 


(*) , 


or *t he 


ccr.trcl 


vector (u) . 




The 


w 


vector 


has zero 


mean value. 


and a covariance 


matrix 


IQb 


wher a: 









E (w) = 0 



( 2 . 4 ) 



and 



id = e[w t : 



( 2 . 5 ) 



The v vector has zero mean value and a covariance matrix 
[R], where: 



E (v) =0 



( 2 . 6 ) 



and 



[B] = E[vv t ] 



( 2 . 7 ) 



The quadratic performance (or cost) index for the linear 
quadradic regulator is the expected value cf: 



2 y"{y T 



J = V 2 / {y T [ A ]y + u r [B]u}dt 



( 2 . 8 ) 



in the statistical steady- state, where [A] represents an 
output cost matrix (a weighting on the output variables); 
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ana [ E ] is a control cost (or control weighting cos 
matrix. 

If full state weighting- is desired, [H] is r = 
by the identity matrix [I]. 



t) 



£ U 



C. PSCGEAS OUTPUT 

1 * C pen-Loo t Eicens vst em Cai cula tic ns 

Ins initial portion of CPIS YSX output includes the 
jrcgiair flew control flags set by the user for that partic- 
ular run, the system of equations being modeled, and the 
cpen-lccp eigenvalue and eigenvector calculations of the [?] 
matrix. 

2 • Bea ulato r Synthesis Ca lcul ation s 

In the solution to the optimal regulator prefclea, 
full state variable feedback is assumed where: 

C C ] = ; [B-‘]£C t ]£SJ (2.9) 



and 



u = -[ c]x (2. 10) 

The control cain [ C ] is a matrix of optimal gains 
which minimize the ccst index expressed in equation (2.8). 

fer optimal regulator synthesis problems, program 
output includes the closed-locp eigenvalues and eigenvec- 
tors; the control gain [C]; the closed-locp dynamics matrix 
[f-GC]; anc the steacy-state gain matrix [S], where [S] is 
the steady-state solution to the algebraic Eiccati equation: 

S[i: + [? T ]S - S[G][E"']£G T ]£ + [ H T ][ A ][ H ] = 0 (2.11) 



13 



a 

ccn tir.ucus 



Kalman filter 
time system may 



cr Estimator which 
be written as: 



a 



e r 



x = [F]x + [K](z - [H]x) 



( 2 . 12 ) 



v r 5 r r • 

•'* .. — s. a 

A 

x is the state estimate 
[K] is a matrix of filter gains, 
and tie state covariarcs is described by: 

I lxx T ) = E jxx T ) + [P* (2.13) 

The filter cain matrix [K] of equation (2.12) is 
obtained from the relationship: 



L n = [P][H T ;[R‘‘] (2. 14) 

••here [ E ] is the steady-state solution to the algebraic 
Eiccati equation: 

[f.'E + P[F T ] - P[F T j[R’ l 3CH]P + [G][Q][G T ] = 0 (2.15) 

, . A 

representing the error covariance of tne estimate x. The 
control covariance is the expected value described by: 



E (UU T ) = [C*x[C T ] (2.16) 

For the Kalman filter/optimal estimator synthesis 
problem, CPTSYSX output includes the eigenvalues and eigen- 
vectors of the optimal estimator (Kalman filter) ; the filter 
gains [K*; the error covariance matrix the covariance 
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sa-rix of tie estimate [x]; the state covariance era 
described in equation (2.13), where 



£ X ] = E£ x x T * 



(2. 17) 



and the control covariance matrix £ 0 ] described in equation 

(2.16) . 

4 . £ t a t io n a ry Clcse d-L ccc Calc ul at ion s 

The stationary response of both state and control 
are presented as roct-mean-sguare values of the state and 
control covariance matrices£X] and £U] described in equa- 
tions (2.17) and (2.16) respectively. 



E. SCIOTICN &LGCRITE2 

Cne cf the fundamental techniques necessary to quadratic 
synthesis cf optimal control systems is the steady-state 
solution cf the algebraic Riccati equation. This is a ncn- 
trivial task due to the iterative nature of the solution. 

The steady-state solution by any quadrature method 
necessitates time increment selection no greater than seme 
fraction cf the shortest period of the closed-loop system; 
imposing a significant computer solution time expenditure on 
the user as well as the requiring an extensive amount of 
computer storage capability due to the matrix expansion 
factor involved. Further, the possible necessity cf a 
time-varying solution of these equations for optimal open 
loop ccrtrcl or estimation requires the inversion cf an 
n X n matrix for each time increment where an unsteady solu- 
tion is desired. 

fl powerful and efficient alternate method cf scluticn 
was develcped by 3ryscn and Hall [Sef. 7], based cn eigen- 
vectci decomposition of the eigensystem of the constant 
coefficient EOLER-L AGEANGE equations. 
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Icz the optimal regulator, these equations take the 
form : 



X 

• 


1 


r 

i *' 


b £ 




1 


1 

v 

1 


X 




4 


-f t 






X 






L 








_ _ 



For the optimal falter or estimator, the equations are 
of the form: 



X 



£ Gd EC ^ j d EH 




V 


h t r" 1 h -f t 




X 


- 
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CPTSXSX, and all earlier program versions, use the 
method cf eigenvector decomposition or the ZOLER-LAC-FAUGE 
equations described in [Ref. 7] for quadratic synthesis cf 
contrcl systems. Program calculations are based or. the QE 
algorithm cf Francis, modified by Wilkinson [Ref. 8]. 
Important features cf this method of eigensystem analysis 

I 

are its improved accuracy, and its insensitivity to widely 
separatee eigenvalues. 

A mere detailed description of the QR algorithm and its 
numerical application to eigensystem analysis may be found 
in [Ref. 8 3 • 



I. PROGRAM OVERVIEW 

CPISYSX and its hi subroutines may be divided into 
three basic categories: 

1) setup ard computation sequencing 

2) data input 

3) calculation 

A brief and general description of the program modules 
and subroutines supporting these basic categories concludes 
this section. 
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1 . 



Pro ble m Eescr iut io n /Program Flew Control 

The MAIN program allows interact ivs selection of all 
program flcv. ccr.trcl flags, and is aided by three s air ra- 
tines : EEEEAL, RDINI, ana RDCHAE. A detailea description 

cf these subroutines is provided in Chapter III. Subroutine 
CEECK verifies the consistency of all reguested program 
eptiers. 



Interact ive Eata Input 



Interactive cata input is performed by the 14 REAl 

subroutines. A detailed description cf the operation of 
these subroutines is also included in Chapter III. Internal 
data generation or external data input is provided by 
subreutire SETUP. 



3- Cal cula tion Seg u enc ing 

Subroutine IEEER functions as a second MAIN program. 
It orders tie data ir put/ca Iculation sequences for the type 
cf problem being solved and performs numerous matrix calcu- 
lations. It is from this subroutine that all input and 
calculation sequences are ordered and performed. 

4 . £3 Alg o r it hm Tra nsf o rmation 

Subroutines MINV, BALANC, OHTHE5 , OUTRAN, HQB, KQE2 , 
EAL3AK, CNCEii, and EREXIT perform the major calculation 
sequences cf the "QR" algorithm. 

5» Eiccati Equation Calculations 

Subroutine EGAIN separates the eigenvalues and 
eigenvectors of the Euler-L agrange equations by eigenvector 
decc I pcsiticn. RGAIN and subroutine SCOV calculate the 
steady-state soluticr cf tie Riccati equations fer the 
controller cr estimator problem. Subroutine SCOV computes 
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the covariance matrix solution to the algebraic 
equation . 



.cam 



6 . ]!cc al Ca 1c ui atlcr.s 

Subroutine HCEZ confutes the modal transformations 
required for modal analysis. 

1 Transfer Function Calculations 



Subroutines IF, FOLES, ZEROS, RESID, ACOMr, CCOMP, 
EQR, anc Function SCI perform transfer function computations 
associated with Modal calculation sequences. Subroutine 
FSDCSI computes the power spectral density of the outputs or 
controls cf a controlled system. 



8. lata Out pu t 



object 

format 



Subroutine S SPENT prints all program calculat 
time format. Subroutine MIT PET allows v 
screen viewing cf all interactive matrix data 



icr.s in 
ariatle 
input. 
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2 . 



EIS3GU CONSIDERATIONS 



Curing the development cf CPTSYSX, all program modifica- 

" iO L S 32 C aadi cions were 2 CCUSed pmma my ZCWard inxersC— 

tive user operation. Experience has demonstrated zhat 
interactive computer cc mmun icaticn offers many advantages in 
the research and pretiem solving environment. The opportu- 
nity for flexible and immediate computer communication , as 
well as the ability to select alternate solution methods, 
are significant advantages to the user; advantages which are 
unavailable in a batch processing environment. 

Although previous versions cf OPTS Y3 produced all the 
desired calculations in the hatch environment, the input 
format ana data sequencing and naming conventions were 
confusing tc many users. The user was burdened with the 
rscessizy cf verifying the correctness of input data format 
and pregram flow control flag seizings for each program run, 
in cider tc ensure the desired calculation sequence was 
properly performed. 

These requirements, combined wizh incomplete pregram 
documentation, promoted a lack cf confidence in the results 
and discouraged many potential users from even attempting to 
use CPTSYS. 



E. EECGE2E LANGUAGE 



CPTSYSX is programmed in FORTRAN IV, following the 
conventions of IEM System/360/370 FORTRAN IV language. Very 
few pregram features have heen been incorporated which are 
ret written in ANSI Standard FORTRAN. These subzle differ- 
ences nctwithst andin c , OPTSYSX has been compiled and run 



under both FORTRAN IV (Gl) and FORTRAN H (Extended) 
compilers cr the I £ l J . 5033. Although tide ever 3 1 _ r rcgiam 
length ns ir. excess cf 2800 lines of text, it ns considered 
completely tcrtahle ficm one operating system to another. 

Cr the presumption that most scientific and research 
perscnnsl are familiar with FORTRAN language, program modi- 
fications may be easily undertaken once system operation is 
understood. 



C. GEKZBAI PROGRAM ECBIFIC ATICNS 

All of the previously developed numerical calculation 
sequences cf OPTSYS were retained un-modified in QPTSYSX. 

Those program sequences requiring the input cf diagonal 
ccst cr covariance matrix elements were deleted cr modified 
to improve user flexibility in entering any desired 
weighting matrix, diagonal or otherwise. This modification 
streamlined program operation through elimination cf several 
program flew control flags; reduced a measure of user uncer- 
tainty; anc further decreased the required degree of user 
program familiarity-promoting uninterrupted operation. 

The FEAE subroutines represent variations on the prin- 
ciple cf simple, effective methods of interactive input, 
coupled with err cr-ccrrecti cn/reccver y sequences. 

Subroutines HDCHAB, RDINT, and RD3EAL were developed by 
the author and Cdr. P.D. Sullivan to accomodate varying 
input format requirements; null-string entry protection was 
developed by Cdr. P.E. Sullivan. These program features are 
discussed in greater detail later in this chapter. 
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E- 



ISIIBACUVE EROGFAE DEV ELCEHENT 



1 • Eicqram Flew Ccr ~.ro 1 

Initial program development required an under- 
stand in c c i the vaiicus problem descriptions as well as 
program input and calculation sequence. After careful anal- 
ysis cf £Hef. 2], a basic program flow control diagram ias 
estatlisred. From this flew control diagram, a logical 
branching network was formulated whose path could be deter- 
mined through either binary logic or numerical selection. 

Three basic branching categories were established 
from the various problem description statements: 



Logical {"Yes" cr "No"..} 

Integer {" 1 « , "2 ", et c. . .} 

Real Number {"Input the value of...} 



Frcm the viewpoint cf free-format computer ccmmuni- 
caticr, integer and real number input presented no signifi- 
cant problems. FOETRAN compiler language is written such 
that numerical data input (real number or integer) is 
expected, thereby requiring only an INTEGER or REAL data 
type statement within the program. Once the data type has 
been declared, the desired values may then be input with a 
free-fermat READ (5,*) statement. 

Cne note of caution concerning numerical data input 
in free-format deserves mentioning: Although all FCETBAN 
compilers treat character string input as an illegal data 
type conversion, many will automatically convert the inad- 
vertent character entry to a "zero" and continue execution. 
Erotecticn against inadvertent errors cf this type is 
discussed later in this chapter. 

Logical input ("Yes" cr "No") poses a unique prcblem 
tc preqiammers. The usual method of incorporation is to 
require the user tc convert the logical answer into an 
icxeger i.e., "Yes" = 1, "No" = 2. These integers may then 
be read directly intc the program, determining program flew. 



Although this method may promote program tiro ease, 
it requires the user to adopt at unnatural hart - , 
pattern — one which -increases the possibility of abnormal 
program termination ir the event of i r.advsrtanf tsar error. 

ii mere refired (from a programming language stand- 
point) arc ergonomic (from the user viewpoint) method of 
logical selection involves utilization cf th 2 entire char- 
acter string answer as an input value. This method has beer, 
inccrpcrated intc 0F1SYSX. 

The logical strings ("Yes" and "No) are declared as 
character strings ir. a data type statement within the 
program cr subprogram. A format statement is also included 
in the program or subprogram utilizing the "A"-field to 
specify tie desired character field width. A REWIND state- 
ment is tier, incorporated in the specific program or subpro- 
gram immediately pricr to each logical string input pcir.t. 
This REWIND statement allows the input device (the terminal 
screen in this case) tc read the character string in the 
same manner as fres-fermat data input. The character field 
width fer this mcdificaticn was established at A1, aliewirg 
streamlined operaticr with the user typing either "Y" cr 
"Yes" fer an affirmative reply; "N" or "No" for a negative 
reply . 

2 • General Input S egu e nc i ng Regu i rements 



All data input to OFTSYSX is in matrix or vector 
format. This data input must be correlated in accordance 
with the problem description and then properly sequenced in 
erder fer the program tc perform the desired calculations. 

The original and modified GPTSYS programs [Ref. 2] 
and [Ref. 3] required not only problem description knowledge 
but complete user familiarity with the detailed calculation 
sequerce cf the program. The latter point was considered a 
significant disadvantage. Elimination of this disadvantage 



was 



<r. s r € a w her a interactive programming offered the 
greatest benefits to the potential user; ans it was toward 
this end that the remaining mcdif icaticns of OPTS'iSX war 
directed. 



3 • Inte rac tive Cat a Inp ut 



In its calculation sequences, OPTSYSX requires the 
input cf up to 14 unique matrices cr vectors. C nee the 
previously described program flow control diagram was 
constructed, data entry points for each matrix cr vector 
were established. At each cf the 16 program data entry 
points, the required input matrix or vector was determined, 
fourteen input subroutines ware added to the original 
program in order to accomodate interactive data entry. 

Ibese matrix input subroutines were written such 

that the user is first informed which specific matrix or 

vector is required; then prompted for the individual matrix 

element values. These values are then individually and 

# • 

sequentially entered from the terminal. Once the matrix or 
vector is filled, if is returned to the terminal screen fer 
user verification and correction if necessary. 

Interactive sequential data entry was programmed by 
means cf a two-dimensional EC loop, with a terminal prompt 
cf the matrix name arc element position prior to the element 
value entry. Data element input is via a f ree-f ormat READ 
(5,*) program statement. 

Cnee the matrix data entry sequence is complete, 
that input matrix is returned to the terminal screen in 
variable format for user ease in row identification. With 
an arbitrary data field width cf 12 characters, the maximum 
number of matrix elements available on an 80 column terminal 
screen is six. Provided the matrix column dimension is less 
than six, this restriction presented no programming fermat 
limitations . 



Per these matrices whose cciumn dimension exceeded 
six, elements were progressively written on sure sc-; :-n~ 
terminal lines. Ones the matrix row is filled t h e screen i? 
dcu tie- spaced, and element display is repeated rn the same 
fashion. This method allows the user to view the matrix 
much as he would expect to see it, providing the advantage 
of ease in row and cciumn identification. 

Kit bin CPTSYSX, subroutine iiATPST performs this 
varratle-f cimat print sequence. The print sequence was 
arbitrarily terminated with a matrix size of 16 X 16, 
presuming that users with iarger systems of equations would 
select alternate forms of data input. 



4. Sa v ir.a Interactive Input 



In most control system design problems, ths system 
matrices generally remain relatively unchanged for a desired 
sequerce cf design calculations. 

In order tc relieve the user of the burden of 
repeated system entry in the interactive mode, several addi- 
tional program flow control flags were added, allowing the 
eptier cf saving the entire original system cf matrices for 
subsequent calculations. Separata options for saving each 
system matrix are automatically offered at the end of each 
program run. 

These matrix saving options provide a further advan- 
tage to the user in that the matrices are redisplayed for 
verification pricr tc calculation execution. Users may then 
change individual matrix elemerts, relieved of the burden of 
full system re-entry. 

5 . Use r-De f inea In put Piles 

Although the basic objective of this work is to 
provide the user with a totally interactive method of data 
input, several disadvantages tc the method of individual 
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' r -? c- 



matrix element input are apparent — input cf very 
laiiiCcc is unwieldy and ti me-ccnsuning ; input cf 
raatricss whose elements remain unchanged from ru: 
iner ficiert . 

In order tc provide an increased measure cf user 
flexibility in data input, subroutine SETUP was modified to 
include previsions for matrix data input from a data file on 
the user's disk. The three system matrices, [F],[G], and 
[ GAM M 2 ] may new be input from the user’s disk. Minor 
program modification is required of the user as follows: 

a. F3ICMS fileaef commands must be modified or 
added tc reflect the name and location cf 
each data set. 

b. The RE2E Format statement (or statements) must 
be changed to reflect the proper data format 
of the user's input file. 



6. Internal Data Generation 



As a further measure cf flexibility, the documenta- 
tion withir subroutine SETUP was expanded to include several 
specific examples of internal matrix data generation. The 
three system matrices [F], [G], and [GAMMA ] may generated 
either within user-written twc-dimensicnal DO-locps or by 
direct assignment statements. These methods may be prefer- 
able for the input of very large matrices with few nen-zero 
element s . 

2 specific example cf internal program generation of 
the output equation [H] matrix is included in subroutine 
REABE . This matrix input method may be preferable for the 
entry cf a large output equation matrix with very few non- 
zero elements. 

Cnee these modifications have been made to subrou- 
tine EE l CP or subroutine READE (as desired) , the program 
should he re-compiled and then run in the usual manner. An 



interactive program ficmpt is provided at the beginning of 
CPTSYSX offering the user the option of specifying oh - 
desired method of data input. 

CPI1YSX vine further nodi find to rr.citr ~ th . s~ flit 7 
to input the [Hj matrix (or ether required input matrices) 
from separate data files. Users with rudimentary program- 
ming skills may now modify subroutine REACH (or one cf the 
ether specific BEAD subroutines) in the manner previously 
described for subroutine SETUP or subroutine REACH. 
Detailed examples of the nature and extent cf these modifi- 
cations may be found in Appendix A. 

7 . Data Entry Corr ecti on 

In an effort tc protect users from errors in data 
input, an error correction sequence was incorporated into 
each matrix input subroutine. 

Cnee the entire matrix or vector is displayed on the 
terminal screen the user is prompted with the questic r. f "Do 
you wish tc change tie value cf any matrix element? Type 
•Yes' or 'No'." If the user types "No", program execution 
continues . 

If the user types "Yes", he is then prompted with 
three additional statements specifying the row number of the 
element tc te changed, the cclumn number of the element to 
he chanced, and the value to be inserted into that matrix 
element. After the corrected value is entered, the new 
matrix values are returned to the screen for 
re-verification. 

This correction sequence continues indefinitely 
until the user signifies that no additional changes are 
recessary. Program execution then proceeds normally. 
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0S2S-2RECR PROTECTION F EATORES 



X a n y 






zter act ivi computer programs suffer tnc 

of chroma! program termination {cini out 
recovery i) should the user inadvert entlv make an erroneous 
keyboard entry. Examples cf obese inadvertent errors 
incl u cs--sc try of a keyboard character or character string 
when the program expects a numerical value; entry cf a 
numerical value when the program is expecting a character 
string; entry cf a null string. In order to preclude 

abnormal program termination dus to these types of inadver- 
tent user errors, several program protection features were 
inccrccrated into OP12YSX. 



1 . Cat a T_ype Conve rsio n Errors 



Three subroutines — EDEEAL, RDINT, HDCHAR — were added 
to OPTSYSX in order to ensure that the proper input data 
type is provided tc the program. Subroutine RDRZAI is 
called at any point a real number or zero integer incutcay 
be encountered; subroutine EEINT is called at any point a 
non-zerc integer input is required; subroutine RECHAR is 
called at any point a character string ("Yes'* or "No") input 
is reguired. 

Kithin each cf these subroutines a null string entry 
protection loop was incorporated (allowing one recovery); 
prompting the user for the correct data type input, and 
returning an error message in the event an incorrect data 
type is encountered. 

fiithin subroutine EDINT, improper data type entry 
was further protected by the addition of a three-way IF 
comparison cf entry integer magnitude. This modification 
precludes illegal (but automatic, with some compilers) data 
type conversion errors. 
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These program design features boast the additional 
advantage cf allowing normal program termination at any 
point in the data input pease by merely pressing toe "Enter" 
key twice. 



2. Inc ons is te nt Erggra m Control Flag Erro rs 



Earlier versions of OEISYS [Ref. 2] and [Ref. 3], 
did provide user error messages in the event cf inconsistent 
program flew control flags, but terminated the program. 
This feature was undesirable from the standpoint cf smooth 
interactive program operation, and was improved in OETSYSX. 

Subroutine CEECK was modified tc include RETURN 
statements any time inconsistent program flow control flags 
are encountered. The user is notified cf the type cf error 
encountered; that run termination has occurred; and prompted 
regarding his desire to return to the beginning of the 
program or terminate execution completely. 
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IV. ESCGBAH DSE AND EXAMPLES 

This chapter contains several basic examples c: tn? 

numerous types of problems which, may be solved using 
C2T3YSX. Included with these examples are copies cf each 
recorded terminal session. 

Potential users should examine carefully the example of 
program failure found in Section D. This example clearly 
demonstrates that unstable mcdes or incorrect choice of 
certain design parameters may cause program failure (and 
incorrect output!), even with the highly stable "QS" algo- 
rithm. It also indicates one possible method of correcting 
this type cf failure by merely making a very small change to 
ere cf the design parameters. 

A- CEEN-LCCP EIG23S ISTE3 A B ALYSIS 

Ihe fcllowing open-loop eigenvalue example was taken 
from [Eef. S, p - 66 9 

Examination of the following program output shews epen- 
lccp eigenvalues at - 1 , -2, and -3. Bote that the eigenvec- 
tors cf the left anc right eigenvector matrices (pa. 33) 
correspond in column fashicn tc the oper.-lcop eigenvalues 
calculated immediately above them (pa. 32). 

The full terminal session is recorded below, with user 
input in lower case le tters following each . 



record cn 

E2GIN BECCB DING CF TEEKINAL SESSION 

E; T=C.O 1/0.02 2 1: 49 :20 

filedef C6 term (reefm fa blksize 133 

global txtlib fortmocl mcd2eeh imsldp ncnimsl 

lead cptsysx (start 
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EXECUTION BEGINS... 

CPTSYSX IS A COMPLET ELY INTERACTIVE OPTIMAL SYSTEMS 
EECGEAK. IT til L I SOLVE NUMEROUS CONTROL P503LZES 
FCLLCWING TYPES CE SYSTEMS CONTROL EQUATIONS: 
XDOT = (F)*X + (G)*U + (GAM) * (W + WG ) 



3 <. 
IK 



a E AS 0 S EM ENT E Q UA 1 10 M- - 
Z = (fi) *X + (D) *W 

REGULATOR PERFORMANCE INDEX — 

0 = 1/2 * INTEGRAL <Yr*(A)*Y + Ut* (5) *U) DT 
STATE FEIEEACK GAIN DEFINITION — 

C = - (C) *x 

DC YCU WISH TC CONTINUE? TYPE "YES" OS "NO". 



yes 



— EATA ENTRY — 



AITHCCGH OP IS Y S X IS SPECIFICALLY DESIGNED TO READ 
AIL MATRIX EATA INTERACTIVELY, SEVERAL ALTERNATE 
METHODS ARE A V AIIABLE TC USERS: 

METHOD 1— THE "F","G", AND "GAMMA" MATRICES 

MAY BE READ FECK SEPARATE DATA FILES. 

MET ECC 2 — THE "F", ,, G" / AND "GAMMA" MATRICES MAY EE 

EXPLICIILY DEFINED WITHIN SUBROUTINE "SETUP" 
(NCTE: IN EITHER CASE, THE USER SHOULD OBTAIN A CCPY 
OF THE ERCGBAM LISTING AND EXAMINE 
THE EXAMPLES CONTAINED IN S/R "SETUP".) 

DC YCC WISH TO CONTINUE? TYPE "YES" OR "NO". 



yes 

EC YCC WISH TO INPUT THE "F", "G", AND "GAMMA" 
MATRICES FROM SUBROUTINE "SETUP" IAK THE 
METHOD DESCRIBED CN THE PREVIOUS SCREEN? 
TYPE "YES" OR "NO". 



EC 



GENERAL OPTSYSX OPTIONS: 

OPTION 1 — SYSTEM ANALYSIS WITHOUT 

CEEN-LOC? EIGEN SYSTEM CALCULATIONS. 
OPTION 2 — SYSTEM ANALYSIS WITH OPEN-LOO? 
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EIGEN SYSTEM CALCULATIONS. 

C IT ION 3 -- CEE l'-LO C F EIGEN SYSTEM FOUND 
2ND PROGRAM TERMINATES . 

I "E ” ~ i ATRIX ENTRY FOLLOWS I MH E £1 AT ZL I . ) 
CPTICN 4 -- MCDAL DISTRIBUTION MATRICES COMPUTED 
WITHOUT FILTER OR REGULATOR SYNTHESIS 
CR STEADY-STATE ANALYSIS. 

SELECT AN OPTICS: 1,2,3, CR 4. 

1 

2 

ENTER THE # OF SIAIES (NS) OF THE SYSTEM MATRIX 
" F "-M ATRI X . 



FLAG/EABAMETER SETTINGS FOB THIS RUN ARE AS FOLLOWS: 
IOL IC IR ISS IE ITF 1 ITF2 ITF3 IFDFW IE ICEEUG 
20 00C0 0 0000 

ISST IDSTAE IPSD IYU INCBM IREG NS NC NOB NG 
0 0 0C 0 03000 

ORDER CE SYSTEM = 3 

NUMEER CF CONTROLS = 0 

NUMEFH CF OBSERVATIONS = 0 

NUMEER CF PROCESS NCISE SOURCES = 0 

ENTER THE SYSTEM MATRIX "F" -MATRIX 

DIMENSION = # STATES (NS) X # STATES (NS) 

TEE ELEMENT F( 1, 1) = 

0 

THE ELEMENT F ( 1 , 2) = 

2 

1 

TEE ELEMENT F( 1, 3) = 



3 1 



0 



TEH ELEMENT F ( 2, 1 ) = 



-6 



-11 



-6 



TEE ELEMENT F ( 2, 2) = 



LEE ELEMENT ? { 2, 3)= 



TEE ELEMENT F ( 2, 1) = 



THE ELEMENT F ( 3, 2) = 



TEE ELEMENT F ( 2, 3) = 



THE SYSTEM MATRIX "F"-MATRIX ... 

C.C 1.0C000 0.0 

C.O 0.0 1.00000 

-6.CC00C -II.OCCOO -6.C0000 

DO YCU WISH TO CHANC-E THE VALUE OF ANY MATRIX ELEMENT 
TYPE "YES” CE "NO". 



n 



OPEN LCCP DYNAMICS MATRIX.. F. 

0.0 0.100CE+01 0.0 

O.C 0 .0 0. 1000D + 0 1 

-0. 60 COE + O 1 -0.110CE + 02 -C. 6000D + 01 

OPEN IOCP EIGENVALUES DET(SI-F).. 

-1.00CC0I+OC:-2.00O0CE+00:-3.00000D+00: 
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OPEN ICCP EIGHT EIGENVECTOR MATRIX 

-5. 77 25 C3E- C 1 - 2. 1 82 175D-01 1.C48285D-G1 

5 . 77 35 CED- C 1 U.364 35SD-01 -3 . 1 44355 D-0 1 

-5. - S . / 2 6 / 1 6 D - 0 1 9.4345640— 01 

OPEN ICCP LEFT EIGEN VICTOR MATRIX T-INV.. 

-5. 1 S 6 1 5 2D ♦ CO - 4.330 127D + 00 - 8 . 66 02 54 D- 0 1 
1.3 747730+ Cl 1.833030D+01 4.532576D+00 

9.539382E+C0 1.430 SC9D + 01 4.769696D+00 

ANALYSIS COMPLETE. CO YOU WANT ANOTHER RUN? 

TYPE "YES" 05 "NO", 
to 

CPTSYSX IS NCK TERMINATED 

5; T = C.26/0.89 21: 52 :08 
record elf 

END EECCFDI NG OF TERMINAL SESSION 
E. E EGUIATCR SYNTHESIS 

Its following regulator synthesis example was taken from 
•'Lecture Notes on Advanced Control Systems", by Professor 
C.J. Ccllirs of the Naval Postgraduate School, Monterey, Ca. 
This example involved determination of the optimal regulator 
gains based on an arbitrarily chosen quadratic index; with 
the various system arc cost matrices described below. 

Examination of the extensive program output indicates 
that the optimal regulator gains are: -5.0 and -SQRT(IC.C). 

The algebraic sign of the gains is consistent with the defi- 
nition displayed on the first screen of the program (pa. 
34) . 

The full terminal session is recorded below, with user 
input in lever case lett ers following each "?" . 

ieccrd on 

EEGIS EECOEEING CF T EFHINAL SESSION 
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?; I=C.G1/C.02 13:55:26 

filecsf C6 term (recfc fa biksize 1 33 

global titlib fortmocl nod2eeh imslap nonimsl 

load cctsysx (stari 

EXEC U TIC N BEGINS... 

CPISYSX IS A COMPLETELY I NT ERACTIVE OPTIMAL SYSTEMS CONTROL 
PFCC-FAM. IT HILL SOLVE NUMEROUS CONTROL PROBLEMS CK TEE 
FCLLCWIKG TYPES CE SYSTEMS CCMTSCL EQUATIONS: 

X EOT = (F) *X + (G ) *U + (GAM) * { W+ WO) 

MS AS LEE ME NT EQUATION — 

2 = ( H) =» X + (D) *S + V 

REGULATOR PERFORMANCE INDEX — 

J = 1/2 * INTEGRAL (Yt*(A)*Y + Ut*(B)*U)DT 
STATE FEEDBACK GAIN DEFINITION — 

U = -(C) *x 

DC YCU WISH TC CONTINUE? TYPE "YES" OR "NO". 



£ IT ECU GH OPTSYSX IS SPECIFICALLY DESIGNED TO READ 
AIL MATRIX LATA INTERACTIVELY, SEVERAL ALTERNATE 
EFTECES ARE AVAILABLE TO USERS: 

KETHCE 1 — THE "F»,»G", AND "GAMMA" MATRICES 

MAY EE READ FROM SEPARATE DATA FILES. 
METHOD 2 — THE "F","G", AND "GAMMA" MATRICES MAY BE 

EXPLICITLY DEFINED WITHIN SUBROUTINE "SETUP". 
( NOTE : IN EITHEF CASE, TEE USER SHOULD C3TAIN A CCPY 
OF THE PROGRAM LISTING AND EXAMINE 
THE EXAMPLES CONTAINED IN S/R "SETUP".) 

DO YOU WISH TC CONTINUE? TYPE "YES" OR "NO". 



y 

EC YOU WISH TC INPUT THE "F", "G", AND "GAMMA" 
MATRICES F f CM SUBROUTINE "SETUP" IAW THE 
METHOD DESCRIBED ON THE PREVIOUS SCREEN? 

TYPE "YES" OR "NO". 



V. 
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G E NERAL CFTSYSX OPTIONS ; 

CHIOS 1 -- SYSTEM ANALYSIS WIIKCUI 

OPEN- ICC P EIGEN SYSTEM CALCULATIONS. 
OPTION 2 -- SYSTEM ANALYSIS WITH OPEN-LOOP 
EIGEN SYSTEM CALCULATIONS. 

CITION 3 -- OPEN- LOOP EIGEN SYSTE M FOUND 
AND P BCG 5 AM TEH MI NAT ES. 

("F" -MATRIX ENTRY FOLLOWS IMMEDIATELY . ) 
CITION 4 -- MODAL DISTRIBUTION MATRICES COMPUTE! 

WITHOUT FILTER. OR REGULATOR SYNTHESIS 



03 STEADY-STATE ANALYSIS. 
SELECT AN OPTION: 1,2,3, OR 4. 



rc 



DC YCU DESISE RMS VALUES CF STATE AND CONTROL? 
TYPE "YES" CR "NO". 

CPTSYSX LQE/CLASSICAL OPTIONS: 



OPTION 


1 


-- CPU MAI FILTER 


AND/OR 


REGULATOR 






SYNTHESIS WITH 


no ext: 


ERNAL "C" 






MATRIX INPUT. 




1 


CPTION 


2 


-- OPTIMAL FILTER 


AND/OR 


REGULATOR 






SYNTHESIS WITH 


EXTERNAL "C" 






MATRIX INPUT. 






CITION 


3 


-- OPTIMAL FILTER 


AND/OR 


REGULATOR 






SYNTHESIS WITH 


EXTERNAL "K" 






MATRIX INPUT. 






CPTION 


4 


-- OPTIMAL FILTER 


AND/OR 


REGULATOR 






SYNTHESIS WITH 


EXTERN 


AL "C" AND 






MATRIX INPUT. 






SELECT 


AN 


CETION: 1, 2, 3, 


OR 4. 





1 



£0 YCD NISH TO DETERMINE THE STEADY-STATE RSSPONSr 
FCE A CONSTANT DISTURBANCE? 
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TYPE "YES" CP. "NO" 



EC 



EC 



YCU WISH 10 DETERMINE THE HOD AL 
AND GAIN MATFICES? 

TYPE "YES" CP "NO". 



DIS 



IBUTICN 



OPEN-LOOP TRANSFER FUNCTION OPTIONS: 

CFIICN 1 — NC OPEN-ICCE TRANSFER FUNCTIONS COMPUTE 

CETICN 2 -- PCLES, RESIDUES, AND ZEROS COMPUTED. 
CHICS 3 — OSH POLES AND ZEROS COMPUTED. 

CETICN 4 — 0 MY POLES AND RESIDUES COMPUTED. 

SELECT AN CETICN: 1, 2, 3, OR 4. 



1 



1 



NOISE TRANSFER FUNCTION OPTIONS: 

OPTION 1 -- NC NOISE TFANSFSR FUNCTIONS COMPUTED. 
C FT ION 2 -- PCLES, RESIDUES, AND ZEROS COMPUTED. 
CEIION 3 -- CNIY POLES AND ZEROS COMPUTED. 

CET ION 4 -- CNLY POLES AND RESIDUES COMPUTED. 
SELECT AN CETION: 1, 2, 3, OR 4. 



COMPENSATOR TRANSFER FUNCTION OPTIONS: 

OPTION 1 — NC COMP. TFANSFER FUNCTIONS COMPUTED. 
OPTION 2 -- ECLES, RESIDUES, AND ZEROS COMPUTED. 
CEIION 3 — CNIY POLES AND ZEROS COMPUTED. 

CETION 4 -- CNLY POLES AND RESIDUES COMPUTED. 

(NOTE: A COMPENSATOR TRANSFER FUNCTION CAN EF 
COMPUTED CNLY IF BOTH A REGULATOfiAND 
FI ITER ARE SYNTHESIZED AND/CR INPUT.) 
SELECT AN OPTION: 1, 2, 3, OR 4. 



1 



Mil A FEED -FOES AED DISTRIBUTION MATRIX 
(" E " - MATRIX) BE INPUT ? 
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TYPE "YES" OE "NO". 



no 



’THIS CPIICN DETER HI S ZS THE CRITERIA FOB DECIDING V HEN A 
HA3KC7 InE AMEX 2 F IS ZZRO-THI MA3K07 PAR A ME IE B IMDICA.ES 
THE CBDEE CF THE NUMERATOR POLYNOMIAL OF EACH TRANSFER 
FUNCTION. 

All "»« 2EH0S OF THIS P CLYNCMIAL ABE PRINTED OUT AND 
THIS TEST TEILS ECS MANY EXTRA ROOTS EXIST AT Z = 0. 
IESS THAN 10.0** (-12) IS CONSIDERED ZERO. 

TEE DEFAULT VALUE CF THIS PARAMETER (IE) IS 6. 

IK CTHEE WORDS, IE = 1.CE-6. 

IF YOU DESIRE A DIFFERENT MARKOV CRITERIA, TYPE THE 
INTEGER VALUE. 

IF YCU DESIRE THE DEFAULT VALUE, TYPE "0" (ZERO) 



C 

DO YCU DESIRE TO SYNTHESIZE A STABLE FILTER (OR REGULATOE) 
BY DESIAEIIIZI NG THE ORIGINAL SYSTEM? 

(NCTElWORKS F CF FILTER OR REGULATOR BUT NOT FOR EOTH 
IN IHE SAME RUN.) 

TYPE "YES" CR "NO". 

EO 

DC YCU DESIBE TC PRINT THE EULER-LAGRANGE EIGEN SYSTEM 
PRICE TO DECOMPOSITION ( EOR CHECKING THE PROGRAM)? 
TYPE "YES" CR "NO". 

yes 

POWER SPECTRAL DENSITY (PSD) OPTION 1 : 

C P IIC N 1 — COMEUTE THE PSD OF THE OUTPUTS AND/CR IHE 
CONTBCLS OE TEE CONTROLLED SYSTEM VEEN 
FORCED EY PROCESS AND MEASUREMENT NCISE. 
(NOTE: ECTH A REGULATOR AND A FILTER MUSTEZ 

RESIDENT IN THE PROGRAM TO USE THIS CPTICK. ( 

C E IIC N 2 — SAKE AS OETICN 1 ABOVE BUT ONLY PRINT TEE 
RESIDUES CE EACH TRANSFER FUNCTION 
USEE IN THE PSD COMPUTATION. 
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CETICN 3 -- NOT DESIRED 
SELECT AN OPTION; 1 



9 ~ J 



OR 3 , 



EC YCU DESIRE REGULATOR SYNTHESIS ONLY? 
TYPE "YES" CR "NO". 



yes 



ENTER THE # OF STATES (NS) OF THE SYSTEM MATRIX 
("I"- MATRIX) . 



r-NT ER TEE # OF CONTFOLS { NC) OF THE CONTROL SYSTEM KCEEI 
( "G" - K A TRI X) . 



ENTER THE # OF MEASUREMENTS CR OBSERVATIONS (NO) OF THE 
( " fi" - i! AT BIX) . 



ENTER IfcE * OF PROCESS NOISE SOURCES (NG) OF THE 
("GAMt^" -MATRIX) . 



C 



FI2G/FAFAMETER SETTINGS FOR THIS RUN ARE AS FOLLCNS: 
IOL IC 15 ISS IM ITF1 IT E2 ITF3 IFDFW IE IEEE UG 
10C0C0 0 0 0 01 

ISSI IDSTAE IPSD IYU INCRM IRSG NS NC NOB NG 
0 C OC 0 12120 



ORDER CF SYSTEM = 2 

NUMEER CF CONTROLS = 1 

NUMEER CF OBSERVATIONS = 2 

NUMEER CF EROCESS NCISE SOURCES = 0 

ENIER THE SYSTEE MATRIX "F"-MATRIX 

DIMENSION = # STATES (NS) X # STATES (NS) 
TEE ELEMENT F( 1, 1) = 
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THE ELEMENT ?( 1, 2 ) = 



1 



I E £ ELEMENT F( 2, 1) = 



0 



TEE ELEMENT F ( 2, 2) = 



C 



THE SYSTEM MATRIX "F"- MATRIX ... 

C.O 1.0CC00 

C.C 0.0 

DC YCC KISH TC CHANGE THE VALUE OF ANY MATRIX ELEMENT? 
TYPE "YES" CR "NO". 



CO 



CEE N LCC? 
O.C 
0.0 



DYNAMICS MATRIX 
0.1000E+01 
0 .0 



F 



CEE N LC CP 

0 . 0 : 



EIGENVALUES 

0 . 0 



DET (SI- F) . . 



CEE N LCCP 
1. 0C000CE+00 
0.0 



SIGHT EIGENVECTCE MATRIX T 

- 1 . 0 00 C C CD + 00 
2 .22 0 4460-16 



C EE N LCCP 
1 . CCCCCCD+CO 
0 . 0 



LEET EIGENVECTOR MATRIX 
4 .5 03 6C0D + 1 5 
4.5Q36C0D+15 



T-INV. . 



ENIER THE MEASUEEMENT SCAIING MATRIX »H"-MATRIX . 

DIMENSION = # OBSERVATIONS (NO) X # STATES (NS) 
TEE ELEMENT H( 1, 1) = 



1 
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THE ELEMENT H ( 1 , 2) = 



r 

V- 

THE ELEMENT H( 2 ., 11 = 

m 

o 

IHE ELEMENT H ( 2 , 2) = 



TEE MEASUREMENT SCALING MATRIX "H"-MATRIX ... 
1.CCCCC 0.0 

0.0 1.0CC00 

EC YOU WISH TO CHANGE THE VALUE OF ANY MATRIX ELEMENT? 
TYPE "YES" OR "NO". 

EO 

MEASUREMENT SCALING MATRIX H. . 

0. 10 COE + Q 1 0.0 

O.C 0.100CE+01 

• 

ENTER THE OUTPUT MEASUREMENT COST MATRIX "A" -MATRIX . 
El ME N SION = * OBSERVATIONS (NO) X # OBSERVATIONS (NO) 
IE E ELEMENT A ( 1, 1) = 

25 

THE ELEMENT A ( 1 , 2)= 

C 

TEE ELEMENT A ( 2, 1) = 

0 

m 

o 

THE ELEMENT A ( 2, 2)= 

C 

TEE OUTPUT MEASUREMENT COST MATRIX "A"- MATRIX ... 
25.CCC0C 0.0 



4 C 



C.O 0.0 

CC YCC KISH TO CHANGE THE VALUE OF ANY MATRIX 



TYPE "YES’* CR "NO". 

OUTPUT COST MATRIX.... A 

C. 2 5 C0E + C2 0.0 

0.0 0.0 



TER THE C CM 13 Cl 



DIMENSION = 4 



1ST FI50TICN MATRIX "G"- HAT FIX . 
STATES (NS) X # CONTROLS (NC) 



TEE ELEMENT G ( 1, 1) = 



0 

THE ELEMENT G ( 2, 1)= 



TEE CONTROL DISTRIBUTION MATRIX "G "-MATRIX ... 

C.O 

1. CCOOG 

DC YCt KISH TO CHANGE THE VALUE OF ANY MATRIX ELEMENT? 
TYPE "YES" CR "NO". ( 

ENTER IHE CONTROL COST WEIGHTING MATRIX "B"-MAIRIX 
DIMENSION = # CONTE CIS (NC) X # CONTROLS (NC) 
IEE ELEMENT B( 1, 1) = 



1 

THE CO NT EC I COST MATRIX E... 

1 . 00000 

DO YCC KISH TO CHANGE THE VALUE Or ANY MATRIX ELEMENT? 
TYPE "YES" CR "NO". 

CO 

TEE CONTROL DISTRIBUTION MATRIX G. . 

0. C 

0. 1CC0E + 0 1 



4 1 



TEE CONTROL COST MATRIX 



B 



0- 1CCCE*C1 



EUIEE-LAGRANGE SYSTEM MATRIX... 



o 

« 

o 


1 .0 00 CCOD + OO 


o 

• 

o 


O 

• 

o 


o 

• 

o 


o 

• 

o 


o 

• 

o 


-1 . C0000CC + Q0 


2.5CCQCCD+C1 


o 

* 

o 


O 

• 

o 


O 

• 

o 


o 

• 

o 


0.0 


- 1 . C00000D+00 


o 

• 

o 



EIGENVALUES AND EIGENVECTORS Or THE 2N X 2N 
EUIER-L AGSAN GE SYSTEM ARIES HQR 2 

- 1 . 5 6 1 1 3 9D+ CO 1.581139D+00 
-1.58 113SE+CO-1. 581 13SD + Q0 
1.58113 9D + CO 1 .581 13SD + 00 
1. 58113SE + C0- 1. 581 13SD+00 



-7. 4 3 C4 4 3D- C2 
4. 1367551- C3 
-1. 154172D+C0 
-3. 5844C6E-C1 



7.16861 2D-02 
- 2 .30 8 345D-0 1 
-2 .0 68 377D-Q2 
-3.7152 2 2D-0 1 



-1 .824925D-01 
-5. 08 24 59D-C2 
2 .63 1337D+00 
-7.5 174 12D-01 



-1. 503482D-0 1 
-5. 262675D-01 
-2. 54 1 22 S D-0 1 
9. 124S27D-01 



EIGENSYSTEM CF OPTIMAL REGULATOR 



EIGENVECTORS FRCM RGAIN PRIOR TO CNORM 
-1. 45SS25D-C1 -2 .6 163 13D-03 
2.34S712E-01 -2 . 266 977D-0 1 

C-LCCE CPTIMAL REG. E- VALUES. .. DET (SI-F+G*C) . . 
-1. 58114E + 00, 1.581 1 4E + 00: 

C-LCCE EIGHT EIGENVECTOR MATRIX H 

-3. 16227 6D-C1 - 3 . 1 6 2 278D-0 1 
1.0C00CCE+C0 0.0 

CCKTECI EIGENVECTOR MATRIX C*M. . 

-1. 531139 E+ 00 1. 581 1 39D + 00 



4 2 



0.0 



C-LCCr CPT . EEG. LEFT E-VECTOR MA TSIX . . H- INV . . 
1 .0 00 C COD+-00 
-2. 16 227SD+C0 - 1 .0 00 C C CD + 00 



IH£ OPTIMAL FELIBACF GalN CONTROL MATRIX. 
-5.0CCGI+G0 -2.1622E+00 






THE CLOSED LCOF DYNAMICS MATRIX F-G*C. . 

0.0 1.0 00 CC CD+00 

-5.0CCCCCD+CO -3.16 2 278 D + 00 

A N AIYS IS COMPLETE. DO YOU WANT ANOTHER F. UN? 
TYPE "YES' 1 OR "NO”, 
no 

OPTSYSX IS NON TERMINATED 

F ; T=C. 42/1.85 14:03:C 3 
record cff 

END FECCFEING OF TERMINAL SESSION 



C. FILTER SYNTHESIS 

The following Kaliran filter synthesis example was taken 
from "Lectcre Notes on Advanced Control Systems", by 
Professor C.J. Collins of the Naval Postgraduate School, 
Monte re y ,Ca . 

This example involved determination of the optimal 
filter gains of an arbitrary system; modeled nearly identi- 
cally to the previous regulator problem. 

In its present configuration, OPTSYSX program sequencing 
reguires the design of an optimal regulator, prior to 
performing any optimal estimator synthesis. In order to 
comply with built-in program sequencing conventions, and 
circumvent program difficulties which may not be specified 
in the particular system model, optimal filter synthesis may 
he accomplished hy entering the identity matrix [I] in those 



program input sequences reguiring the entry of on curgu - 
cost (weighting) matrix. Although tn<. optimal ..jjQlapc. 
calculations nay differ free those expected, the optimal 
e sti mator calculations will be correct for the system medsi. 

Examination cf the extensive program output .indicates 
that the optimal filter gains are: -5.0, and -SQRT(2.0). 

Its full terminal session is recorded below, with user 
input in lower case le tters following each "?" . 



record cr. 

BEGIN EECOBEING CF TEEMINAL SESSION 

E; T = C. 0 1/0 .02 2 1:49:20 

fileaef C6 term (reefm fa blksize 133 

global txtlib fortmod2 mcd2eeh imsidp ncnimsl 

lead cptsysx (start 

EXEC C TICK BEGINS... 

OPTS YSX IS A COMPLETELY INTERACTIVE OPTIMAL SYSTEMS CC MB CL 
PHCGBAM. IT Hill SOLVE NUMEROUS CONTROL PROBLEMS ON THE 
FOLLOWING TYPES CE SYSTEMS CONTROL EQUATIONS: 

XDOT = (E ) *X + (G) *U + (GAM) * (W+W 0) 

MEASUREMENT EQUATION — 

Z = ( H ) *X + (D) *W + V 

REGULATOR PERFORMANCE INDEX — 

J = 1/2 * INTEGBAL (Yt*(A)*Y + Ut * ( E) *U) ET 
ST AT F FEEDEACK GAIN DEFINITION — 

U = -(C) *x 

EO YOU WISE TO CONTINUE? TYPE "YES" OR "NO". 

yes 

— EAT A ENTRY — 

AIIHCOGH OPTS Y SX IS SPECIFICALLY DESIGNEE TO RFAC 
ALL MATRIX EATA INTERACTIVELY, SEVERAL ALTERNATE 
METHODS ARE A V SUABLE TC USEBS: 

METHOD 1 — THE "F","G", AND "GAMMA" MATRICES 
MAY BE BEAD FRCM SEPARATE DATA FILES. 

METHOD 2 — THE "F","G", AND "GAMMA" MATRICES MAY EE 
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E XF II Cm Y DEFINED WITHIN SUBROUTINE "SETUP” 
(NOTE: IN E IT FIS CASE# THE USER SHOULD OBTAIN A 00 



OF 


THE PROGRAM LISTING 


AND 


EXAMINE 


THE 


EXAMPLES CONTAINED 


T 'I 


S/E 5 


• r> 7 r *1 - <) \ 

. X vj ^ 0 f 


DC YOU WISH 


TC CONTINUE? TYPE 


"YZS 


” cs 


"NO”. 



i'L 



yes 

DO IOC WISH TO INPUT THE ”F”, ”G”, AND “ GAiu-i A" 
MATRICES FECK SUBROUTINE "SETUP” I A w THE 
METHOD DESCRIBED C N THE PREVIOUS SCREEN? 

TYPE "YES" OR "NO”. 

EO 

GENERAL CP1SYSX OPTIONS: 

CPI ION 1 -- SYSTEM ANALYSIS WITHOUT 

CPEN-LO CP EIGEN SYSTEM CALCULATIONS. 
OPTION 2 — SYSTEM ANAIYSIS WITH OPEN-LOOP 
EIGENSY STEM CALCULATIONS. 

OPTION 3 -- CPEN-LO C P EIGENSYSTEM FOUND 
AND PROGRAM TERMINATES. 

(”F”-MATBIX ENTRY FOLLOWS IMMEDIATELY.) 
OPTION 4 -- MODAL DISTRIBUTION MATRICES COMPUTED 
VIIHOUT FI IT ER OR REGULATOR SYNTHESIS 
CB STEADY-STATE ANALYSIS. 

SELECT AN OPTION: 1,2,3# OR 4. 

•? 

1 

£0 YCU DESIRE RMS VALUES OF STATE AND CONTROL? 

TYPE "YES” Cfi "NO”. 



no 



CPTSYSX LQ E/CLASSICAL OPTIONS: 

OPTION 1 -- OPTIMAL FILTER AND/OR REGULATOR 

SYNTHESIS WITH NO EXTERNAL ”C" OR "K” 
MATRIX INPUT. 

OPTION 2 — OPTIMAL FILTER AND/OR REGULATOR 
SYNTHESIS WITH EXTERNAL ”C” 

MATRIX INPUT. 
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CFTIOt! 3 



OPTION 



OPTIMAL FILTER 
SYNTHESIS WITH 
MATRIX INPUT. 
OPTIMAL FILTER 
SYNTHESIS WITH 
MATRIX INPUT. 



AND/OR REGULATOR 
EXTERNAL " K " 

AND/OR REGULATOR 
EXTERNAL "C" AND "K" 



SELECT AN OPTION: 1 , 2, 3, OR 4. 



DC 5CC WISH TO DETERMINE THE STEADY-STATE RESPONSE 
ECR A CONSTANT DISTURBANCE? 

TYPE "YES" CR "NO". 

DO YOU KISH TO DETERMINE THE MODAL DISTRIBUTION 
AND GAIN MATRICES? 

TYPE "YES" CR "NO". 

CPEN-ICCP T FA NS F ER FUNCTICN OPTIONS: 

CFIICN 1 — NO OPEN-LOOP TRANSFER FUNCTIONS CCMFUIEC 
CP1ICN 2 -- POIES , RESIDUES, AND ZEROS COMPUTED. 
CPTICN 3 — ONLY POLES AND ZEROS COMPUTED. 

CFIICN 4 -- CNIY PCLES AND RESIDUES COMPUTED. 

SELECT AN OPTION: 1, 2, 3, OR 4. 



NCISE TRANSFER FUNCTION CFTICNS: 

CETION 1 -- NO NOISE TRANSFER FUNCTIONS COMPUTED 
OPTION 2 — POLES, RESIDUES, AND ZEROS COMPUTED. 
OPTION 3 — ONLY POLES AND ZEROS COMPUTED. 

OPTION 4 -- ONLY ECLES AND RESIDUES COMPUTED. 
SELECT AN OPTION : 1, 2 . 3, OR 4. 



COMPENSATOR TRANSEER FUNCIICN OPTIONS: 

OPTION 1 — NO COMP. TRANSFER FUNCTIONS COMPUTED 



OPTION 2 
OPTION 3 
OPTION 4 
(NOTE 



SELECT 



-- POLES, EES IDUES , AND ZEROS COMPUTED. 

— OKIE POLES AND ZEROS COMPUTED. 

— 0 MY POLES AND RESIDUES COMPUTED. 

: A COMPENSATOR TRANSFER FUNCTION OAR .. 0 
COMPUTED ONLY IF BOTH A REGULATOR 
AND FILTER ARE SYNTHESIZED 
ANC/CE INPUT.) 

AN OPTION: 1, 2, 3, OR 4. 



Mil A FEED -? OR fc ARD DISTRIBUTION MATRIX 
"D" - MATRIX EE INPUT ? 



TYPE "YES” OR "fC". 



r.o 

THIS CETICN DETERMINES THE CRITERIA FOR DECIDING WHEN 
A MARKOV PARAMETER IS ZSRC-1HE MARKOV PARAMETER 
INDICATES THE CRD E F OF THE NUMERATOR POLYNOMIAL OF EACH 
TRANSFER FUNCTION. 

• » 

ALL "N" ZEROS CF THIS POLYNOMIAL ARE PRINTED OUT AND 
THIS I EST TELLS HOW MANY EXTRA ROOTS EXIST AT Z = 0. 

LESS THAN 10.0** (-IE) IS CONSIDERED ZERO. 

THE DEFAULT VALUE CF THIS PARAMETER (IE) IS 6. 

IN CTHER WORDS, IE = 1.0E-6. 

IF YOU DESIRE A DIEEERENT MARKOV CRITERIA, 

TYPE TEE INTEGER VALUE. 

IF YCU DESIRE THE DEFAULT VAIUE, TYPE "0” (ZERO) 

•j 

C 

DO YCU DESIRE TO SYNTHESIZE A STABLE FILTER (OR REGULATOR) 
BY DESTABILIZING THE ORIGINAL SYSTEM? 

( NOT E : WORKS EOR EILTER OR REGULATOR BUT NOT FOR ECTH 
I K THE SAME HUN. ) 

TYPE "YES" CR "NO". 
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1*1 



EC YOU DESIFS TC PRINT THE EULER -LAGRANGE EIGEN SYS TEN 



ISICB TC DECC MP03IT ION (?0 3 CHECKING THE PROGS A t) 
TYPE ."YES" CR "NO". 



PCSER SPECTRAL DENSITY (PSD) OPTION 1 : 

CPT1CN 1 — COMPUTE THE PSD OE THE OUTPUTS AND/OR 
THE CONTROLS OF THE CONTROLLED SYSTEM 
WHEN FORCED E Y PROCESS AND MEASUREMENT 
NOISE. (NOTE: ECIH A REGULATOR AND A 
FILTER MUST EE RESIDENT IN THE PROGRAM 
TO DSE THIS CFIION.) 

CP TIC N 2 — SAME AS OPTICS 1 ABOVE BUT ONLY PRINT THE 
RESIDUES OF EACH TRANSFER FUNCTION 
USEE IN THE PSD COMPUTATION. 

CPTJCN 3 — NOT DESIRED. 

SELECT AN CPTICN: 1, 2, OR 3. 



DC YCU DESIRE REGULATOR SYNTHESIS ONLY? 

TYPE "YES" CR "NO". 

ENTER THE # OF STATES (NS) OF THE SYSTEM MATRIX 
"F"-M ATRIX . 



ENIIR I HE * OF CONTROLS (NC) OF THE SYSTEM MODEL 
"G"-M ATRIX . 

•j 

1 

ENTER THE # OF MEASUREMENTS OR OBSERVATIONS (NO) 
"H"- MATH IX . 

7 

1 

ENIIR THE # OF PROCESS NOISE SOURCES (NG) 

"GAMS A "-MATRIX . 



U8 



c 



FLAG/PAFAMETER S 
1QL 3 C IR ISS IS 
0 0 C 0 

ISST IDSTAE IP SD 
0 C 0 



SITINGS FCF THIS HUM 
I IF 1 IIF2 ITF3 
C 0 0 0 

IYU IN CRH IREG NS 
CO 0 2 



ARE AS FOLIC 7 
IFDFW IS I LI 
0 0 
NC NOB NG 

1 1 0 



03C IF OF SYSTEM = 2 

NUMBER CF CCMTHC1S = 1 

NUMEE5 CF OBSERVATION S = 1 

NUMBER CF PROCESS NOISE SOURCES = 0 



ENTER THE SYSTEM MATRIX "F"-MATRIX 

DIMENSION = # STATES (NS) X if STATES (NS) 
TEE ELEMENT F( 1, 1) = 



0 

IHE ELEMENT F ( 1,2) = 



1 

TEE ELEMENT F( 2, 1) = 

0 

THE ELEMENT F ( 2, 2)= 



C 



TEE SYSTEM MATRIX "F'»-MATRIX ... 



C.O 1.0C000 

C.O 0.0 

DO YCL WISH TO CHANGE THE VALUE OF ANY MATRIX ELEMENT 



TYPE "YES" CF. "NO' 1 . 



EO 



CPEN LC CP DYNAMICS MATRIX 
0.0 0.100CE+01 



F 



as 



O.C 0.0 

ZNIE3 IKE MEASUREMENT SCALING MATRIX n H"-MAT3IX . 

DIMENSION = « OBSERVATIONS (NO) X S STATES (NS) 

TEE ELEMENT H ( 1, 1) = 

1 

TEE ELEMENT H { 1 , 2) = 

✓ 

C 

TEE MEASUREMENT SCALING MATRIX "H" -MATRIX ... 
1.CC0CC 0.0 

EC YCU NISH TO CHANGE THE VALUE OF ANY MATRIX ELEMENT? 
TYPE "YES” CS " NO", 
no 

MEAS CEE KENT SCALING MATRIX H. . 

0. 1CC0E + 0 1 0 .0 

ENTER IHE OUTPUT MEASUREMENT COST MATRIX "A" -MATRIX . 
DIMENSION = if CESERVATICNS (NO) X # OBSERVATIONS (NO) 
TEE ELEMENT A( 1, 1) = 



THE OUTPUT MEASUREMENT CCST MATRIX " A"-MATRIX ... 

1. ccooc 

DC YCU WISH TO CHANGE THE VALUE OF ANY MATRIX ELEMENT? 
TYPE "YES" CR "NO”, 
no 

C C IP Cl COSI MATRIX A.. 

0. 10C0D+0 1 

ENTER THE CCNTBCI DISTRIBUTION MATRIX »G"-MATRIX . 

DIMENSION = # STATES (NS) X # CONTROLS (NC) 

TEE ELEMENT G< 1, 1) = 



0 



5 C 



0 



TEE ELEMENT G( 2, 1}= 



THE CCNIROI DISTRIBUTION MATRIX "G "-M AT R IX ... 
C..C 
0.0 

DO ICC WISH TC CHANGE THE VALUE OF ANY M AT 3 IX ELEMENT? 
TYPE "YES' 1 CR "NO". 

yes 

ENTER IKS SCW ND22E3 OF I BE ELEMENT TO 3E CHANGED. 



z 

ENTER THE COLUMN NUMBER CF THE ELEMENT TC BE CHANGED. 



1 

THE ELEMENT G( 2, 1) = 



1 

THE CONTROL DISTRIBUTION MATRIX "G "-MATRIX ... 

C.O 

1. ccooo 

DC ICC WISH TC CHANGE THE VALUE OF ANY MATRIX ELEMENT? 
TYPE "YES" CR "NO", 
no 

ENTER THE CCNTRCL COST WEIGHTING MATRIX "B"-MATRIX 
DIMENSION = t CON 1 RC IS NC X # CONTROLS NC 
THE ELEMENT B ( 1 , 1 ) = 

1 



TEE CCNTRCI COST MATEIX B... 

1. ccooc 

DC YOU WISH TO CHANGE THE VALUE OF ANY MATRIX ELEMENT? 
TYPE "YES" CR "NO", 
no 
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TEE CONTROL DISTRIBUTION MATRIX 



G 



0. C 

0. 10C0I+0 1 

THE CONTROL COST 2 AT SIX 15. . 

0. 1CCGE+G 1 

EIGENS Y STE M CF CETIMAL REGULATOR 

C-LCCP OPTIMAL 3EC-. E- V ALU E S . . . DE T ( SI- F+G*C) . . 

-7 . C 7 1 ( 7E- C 1 , 7.C7107E-01: 

C-ICCP RIGHT EIGENVECTOR MATRIX M 

-7. 07 1066E-01 -7.071 C6SD-01 
1.CCC0CCD+C0 0.0 

CONTROL EIGENVECTOR MATRIX C*M.. 

-7 .C71066D-01 7.C71068D-01 

C-LCCP OPT. REG. LEFT E-VECTOR MATRIX. . M- IN V . . 

0.0 1 .0 00 CC0D+00 

-i. 4 142 i4E+co - 1 .ooocccd+oo 

THE CFTIMAI FEEE3ACX GAIN CONTROL MATRIX C=3INV*G1*S . 

-1 . OC CCE+OQ -1.414 2E + 00 

THE CIO SED L CO P DYNAMICS MATRIX F-G*C. . 

0.0 1. OOOCCCD+OO 

-1 . OCCCCCD+CO -1 .41 4214D+00 

ANAIYSIS COMPLETE. EO YOU WANT ANOTHER RUN? 

TYPE "EES" OR "NC". 
no 

..OPTSYSX IS NOW TERMINATED 

F; T = C. 54/2.84 1 5: 36:4S 
ieccid off 

END EECCREING OF TERMINAL SESSION 



E. E3AME1E CF PBOGEAE FAIL DEE 

The following pathological sxample cf program tailors 
during regulator synthesis was oaken fro a tne Journal or 
Guidance and Control, Vcl.3, Nc.2, pp. 190-192, March-Acril 
1980. 

Ir. ohis example, the choice of the quadratic index value 
was the factor promoting program instability; leading to 
eventual program failure in subroutine HQR2 . The calculated 
regulator gains cf -5.1, and -3.1 (pa. 63) are not ccrrecti 

With a ’slight* iicdifi cation of the ccst matrix from a 
previous value cf 4.0000 to a new value of 4.0001, the 
program was run a second tine. Failure did net occur cn the 
second run, anc the new calculations (pa. 68) indicate 
filter gains of -2.0, and a "small" residue of 3. 19D-14 
(essentially zero). These are the correct values. 

This example pcints cut one possible method cf 
ccrrectirg certain pregram failure modes, shculd they cccur 
during execution. 

Tie full terminal session is recorded below, with user 
input in lower case le tter s fc Hewing each "?" . 

Following the program failure example, that pertien cf 
the repeated terminal output was delated up to the point 
where program execution of the second run begins. 

record cn 

E2GIN EZCCEEING CF TEfKINAL SESSION 

5; T = 0. 0 1/0.02 2 1: 49:30 

filedef C6 term (reefm fa blksize 133 

global txtlib fortmoc2 mcd2eeh imsldp ncnimsl 

lead cptsysx (start 

EXECUTION EEGINS . . . 

OFTSYSX IS A COMPLETELY INTERACTIVE OPTIMAL SYSTEMS CCNTECL 

PEOGEAtt. IT WILL SOLVE NUMEROUS CONTROL PROBLEMS ON TEE 
FOLLOWING TYPES OF SYSTEMS CONTROL EQUATIONS: 
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XDOT 



= (F)*X + (G)*U + (GAM) * (« + W0) 

MEASUREMENT EQUATION — 

Z = (H ) *X + (D) *W + V 

iiriGU I AxO H PERFORMANCE INDEX - — 

J = 1/2 * INTEGRAL (Y * (A) *Y + U *{B)*U)£T 
STATE FEEDBACK GAIN DEFINITION— 

U = - (C) *x 

DO YOU WISH TO CONTINUE? TYPE "YES" OR "NO". 



ALTECUGH OPTSYSX IS SPECIFICALLY DESIGNED TO HEAD 
All MATRIX DATA INTERACTIVELY, SEVERAL ALTERNATE 
KFTECES ARE AVAILABLE TO USERS: 

MEIHCD 1 — THE "F","G", AND "GAMMA" MATRICES 
MAY BE EEAD FROM SEPARATE DATA FILES. 

METHOD 2 — THE "F'V'G", AND "GAMMA" MATRICES MAY EE 
EXPLICITLY DEFINED WITHIN SUBROUTINE "SETUP'. 
(NOTE: IN EITHER CASE, THE USER SHOULD OBTAIN A COPY 
OF THE EROGRAM LISTING AND EXAMINE 
THE EXAMPLES CONTAINED IN S/R "SETUF". ) 

EC YOU WISH TO CONTINUE? TYPE "YES" OR "NO". 



yes 

EC YOU WISH TO INPUT THE "F", "G", AND "GAMMA" 
MATRICES FFCM SUBROUTINE "SETUP" IAW THE 
METHOD DESCRIBED ON THE PREVIOUS SCREEN? 
TYPE "YES" OB "NO". 



no 



GENERAL OETSYSX OPTIONS: 

OPTION 1 — SYSTEM ANALYSIS WITHOUT 

CEEN-LOCF EIGEN SYSTEM CALCULATIONS. 
OPTION 2 — SYSTEM ANALYSIS WITH OPEN-LOOP 
EIGENSY STEM CALCULATIONS. 

OPTION 3 — CPEN-LO CP EIGENSYSTEM FOUND 
AND EROGRAM TERMINATES. 

( "F"~MATRIX ENTRY FOLLOWS IMMEDIATELY.) 
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CPIION 4 — 2CEAL D 1ST RIEUTION MATRICES COMPUTED 
WITHOUT FILTER 03 REGULATOR SYNTHESIS 
Cr STE'A EY-ST ATE ANALYSIS. 

SELECT AN OPTION: 1,2, 3, 03 4. 



yes 



EC 



EC ICE DESIRE RES VALUES CF STATE AND CONTROL? 

TYPE "YES " CR "NO". 

CPTSYSX LQ3/CLA SSICA L OPTIONS: 

OPTION 1 — OPTIMAL FILTER ANE/OR REGULATOR 

SYNTHESIS WITH NO EXTERNAL "C" OR "K" 
MATRIX INPUT. 

CETION 2 -- OPTIMAL FILTER AND/OR REGULATOR 
SYNTHESIS WITH EXTERNAL "C" 

MATRIX INPUT. 

OPTION 3 -- OPTIMAL EILTER ANE/OR REGULATOR 
SYNTHESIS WITH EXTERNAL "K" 

MATRIX INPUT. 

OPTION 4 — OPTIMAL FILTER AND/OR REGULATOR 

SYNTHESIS WITH EXTERNAL "C" ANE "K" 
MATRIX INPUT. 

SELECT AN CPTICN: 1, 2, 3, OR 4. 



EO YOU WISH TO EETERMI NE THE STEADY-STATE RESPONSE 
FOR A CONSTANT DISTURBANCE? 

TYPE "YES" CR "NO". 

EC YOU WISH TO EETERMI NE THE MODAL DISTRIBUTION 
AND GAIN MATRICES? 

TYPE "YES" CR "NO". 



CPFN-LCOP TRANS I FR FUNCTION OPTIONS: 

CPTICN 1 -- NO CPEN-LCCP TRANSFER FUNCTIONS COMPUTEE. 
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CPTICN 2 -- POLES, RESIDUES, AND ZEROS COMPUTED. 
CP1IC.N 3 -- OSH POLES AND ZEROS COMPUTED. 

CPTICN 4 -- ONLY POLES AND RESIDUES COMPUTED. 
SELECT AN OPTION : 1, 2, 3, OR 4. 



NCI S E TRANSFER FUNCTION OPTIONS; 

OPTION 1 — NC NOISE TRANSFER FUNCTIONS COMPUTED. 

OPTION 2 -- POLES, RESIDUES, AND ZEROS COMPUTED. 
OPTION 3 -- ONLY POLES AND ZEROS COMPUTED. 

OPTION 4 -- ONLY POLES AND RESIDUES COMPUTED. 
SELECT AN OPTION; 1, 2, 3, OR 4. 



1 



■? 

1 



COMPENSATOR TRANSFER FUNCTION OPTIONS: 

OPTION 1 — NC COMP. TEANSFER FUNCTIONS COMPUTED. 
OPTION 2 -- POLES, RESIDUES, AND ZEROS COMPUTED. 
OPTION 3 — C MY POLES AND ZEROS COMPUTED. 

OPTION 4 -- ONLY POLES AND RESIDUES COMPUTED. 

(NOTE; A COMPENSATOR TRANSFER FUNCTION CAN E 
COMPUTED ONLY IF BOTH A REGULATOR 
AND FILTER ARE SYNTHESIZED 
AND/OR INPUT.) 

SELECT AN OPTION : 1, 2, 3, OR 4. 



WILL A FEED-FORiiARD DISTRIBUTION MATRIX 
("D" - MATRIX) EE INPUT ? 

TYPE “YES" OR "NO". 



to 

DO YOU EISIBE TO SYNTEESIZE A STABLE FILTER (OR REGUIATCR) 
BY DESTABILIZING THE ORIGINAL SYSTEM? 



( NCT E ; WORKS FOE FILTER OR REGULATOR BUT NOT FCE EOT H 
IN THE SAME RUN.) 
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TYPE "YES" CE "NO" 



DC YOU DE SI 53 IC 5 HI NT THE EULER-LAGRANGE SIGEN5Y: 
PRICE TO DECCNPOSITT.ON (EOS CHECKING THE PROOF. 
TYPE "YES" CS "NO". 



j. _ i. 

v\ i 



r.c 

PCiiEE SPEC THAI DENSITY (PSD) OPTION 1 : 

CPTICN 1 — COMPUTE TEE PSD OF THE OUTPUTS AND/OR THE 

CONTROLS CF THE CONTROLLED SYSTEM NHEN FORCED BY 
PROCESS AND MEASUREMENT NOISE. (NOTE: ECTH A 
REGULATOR AND A FILTER MUST BE RESIDENT IN THE 
PROGRAM TC USE THIS OPTION.) 

CPTICK 2 — SAME AS CETION 1 ABOVE BUT ONLY PRINT THE 
RESIDUES CF EACH T5ANSFER FUNCTION 
USED IN TEE PSD COMPUTATION. 

CPTICN 3 — NOT DESIRED. 

SELECT AN CITION: 1, 2, OR 3. 



DC YOU DESIRE REGULATOR SYNTHESIS ONLY? 

TYPE "YES" CR "NO". 

yes 

ENTER THE # OF STATES (NS) OF THE SYSTEM MATRIX 
("F»-K ATRIX) . 

* 

2 

ENTER THE # OF CCNTROLS (NC) OF THE SYSTEM MODEL 
( »G"- MATRIX) . 

* 

1 

ENTER THE * OF MEASUREMENTS OR OBSERVATIONS (NO) 
("H "-MATRIX) . 



ENTER THE # OF PF0CZ5S NOISE SOURCES (NG) 
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("G AM K A"- M A 161 X ) . 



C 



IOL 

1 

1 SE1 
0 



ELAG/EAFAMETER S 
IC IE ISS IK 
1 C 0 
ID SI A E IP SD 
C 0 



ETTINGS FOB 
ITF1 IT F 2 
0 1 0 
IYU INCBM 
C C 



THIS sail ARE 
ITF3 IF DEW 
0 0 
IREG NS NC 
1 2 1 



AS FCIIC ..S : 
IE IDEEUG 
0 0 
NOB II G 
2 0 



06DEE CE SYSTEM = 2 

NUKE EE CE CONTROLS = 1 

NUKEEE CF CBSER VATIC NS = 2 

NUKEEE CE EROCESS N C IS E SOURCES = 0 



ENTER THE SYSTEM MATRIX ( «F "-MATRIX) 

DIMENSION = # STATES (NS) X * STATES (NS) 
THE ELEMENT F ( 1, 1) = 

1 

C 

TEE ELEMENT F( 1, 2) = 



1 

THE ELEMENT F ( 2, 1) = 

-1 

TEE ELEMENT F( 2, 2) = 

1 

0 

THE SYSTEM MATRIX ("F"-MATRIX) .. . 

0.0 1.0CC00 

-1.CCOOC 0.0 

EC YCU WISH TO CHANGE THE VALUE OF ANY MATRIX ELEMENT? 
TYRE "YES" Cfi "NO". 

TiC 
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CEO LCCF DYNAMICS MATRIX 
O.C 0.100CE+01 

-0. 1CC0E+G 1 0.0 



CEO ICCP 

0 . 0 



EIGEN V AIDES. 

1. C0Q0CE + 00: 



CST(SI-F) ♦ • 



CEE N LCCF 

0 . 0 

1. CCCCCCD+CO 
CEE N ICC? 



RIGHT EIGZNVSCTCF MATRIX I.... 

- 1 .000 CCCD + 00 

0.0 

LEET EIGENVECTOR MATRIX T-INV.. 



0.0 1.000CCCD+00 

-1. 00CCCCE+C0 0.0 



ENTER THE HEASOEEMENT SCALING MATRIX ("H"-M ATRIX) . 

DIMENSION = # OBSERVATIONS (NO) X 4 STATES (NS) 
THE ELEMENT H( 1, 1) = 



C 

TEE ELEMENT H( 1 , 2) = 

*5 i 

m 

o 

IKE ELEMENT H ( 2, 1 ) = 



C 

TEE ELEMENT H( 2, 2) = 



TEE MEASUREMENT SCALING MATRIX ( " H" -M AT RIX ) . . . 

0.0 0.0 

C.O -1.0CC0O 

IC YCU WISH TO CHANGE THE VALUE OF ANY MATRIX ELEMENT? 
1 IPS "YES" CR "NO". 



no 
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MEASUREMENT SCALING HAT BIX 



H 



0.0 0.0 

0. C - 0. 10001*01 

I'CEAI MEASUREMENT SCALING MATRIX H(3AR)*T.. 

C.O 0.0 

- 1 .CCCOOOD+OO o.c 

ENTER THE OUTPUT MEASUREMENT COST MATRIX ( ” A " - M ATS IX ) . 
II K ENSIGN = # OBSERVATIONS (NO) X * OBSERVATIONS (NO) 
THE ELEMENT A( 1, 1) = 



C 

TEE ELEMENT A( 1 , 2) = 

0 

THE ELEMENT A ( 2, 1)= 

✓ 

C 

TEE ELEMENT A( 2, 2) = 

* 

a 

THE OUTPUT MEASUREMENT COST MATRIX (" A"-MATRIX) . . . 

C.O 0.0 

C.O 4.0CCO0 

IC YOU WISH TO CHANGE THE VALUE OF ANY MATRIX ELEMENT? 
TYPE '’YES" CB "NO". 

BO 

CC1PCT COST MATRIX A.. 

0.0 0.0 

0.0 0.40 00E+0 1 

ENTER THE CCNTECL DISTRIBUTION MATRIX (»G"-MATRIX) . 

DIMENSION = # STATES (NS) X S CONTROLS (NC) 

TEE ELEMENT G ( 1 , 1 ) = 
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c 

TEE ELEMENT G( 2, 1) = 

7 

1 

TEE C C ST SC I DISTRIBUTION MATRIX ( "G "-H AT R IX) . . . 

C.O 

1. cccoc 

EC YCU WISH TO CHANGE TEE VALUE OF ANY MATRIX ELEMENT? 

TYPE "YES" CE "MO", 

no 

ENTER THE CCNIRCI COST WEIGHTING MATRIX ( "3" -M ATS I X) 

CIMENSION = # CONTROLS (NC) X # CONTROLS (NC) 

THE ELEMENT B( 1, 1) = 

7 

1 

THE CONTROL COST MATRIX E... 

1.CC00C 

EC YCU WISH TO CHANGE TEE VALUE OF ANY lNAT3IX ELESENI? 

TYPE "YES" CE "NO", 

no 

TEE CONTROL DISTRIBUTION MATRIX G. . 

0.0 

0. 1CC0E+01 

KCEAI CONTROL DISTRIBUTION MATRIX TI*G.. 

1 .CCCOC CD+ 00 

0.0 



THE CONTROL COST MATRIX.. B 

0. 10C0E + O 1 

OPEN LOOP TRANSFER FUNCTIONS... 

TF FCE INPUT NO. 1 AND OUTPUT NO. 1: 



= 0 . 0 
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NC FINITE ZEROS 



TF GAIN 



RESIDUES AT THE PCIES: 

JC1ES RESIDUES 

REAL { A ) I U AG ( E ) 

( 0.0 ) +J ( 1.0C000C) ( 0.0 ) EXP (A*T) *COS (B*T) 

( 0.C ) + J (-1 . 000000) ( 0.0 ) EXP <A*T) *SIN (B*T) 

I? EC B INPUT NO. 1 A ND OUTPUT NO. 2; 

ORDER CE NUMERATOR = 1 TF GAIN = -0.10000+01 



NUMERATOR EIGENVALUES (INCLUDING EXTRANEOUS ZERO VALUES) 
( 0.0 )+J | 0.0 ) 

( 0. 0 ) + J ( 0.0 ) 

RESIDUES AT THE FOIES: 



POLES 

REAL (A) IMAG (6) 

( 0 . 0 ) +0 ( 1 . 000000 ) 
( 0. C ) +J (-1.000000) 



RESIDUES 

( -1.C00000) EXP (A*T) *COS (E*T) 

( 0.0, ) EXP (A*'I) *SIN (E*T) 



E AILUEZ IN HQE2 OK EIGENVALUE NO. 4 



-1. 962366D + C0 
3.46 483 SD- 03 
-4. 4 1504 ID- 15 
5.28 1S45D- 11 



3.4646 1 2D-03 
3 .762 172D-02 
-3.20 8 S43D- 1 3 
- 1 .2 67 8 12D-1 7 



-2.4S9867D+00 
-1 .49 1143D+00 
- 1.962366D+00 
3 .621125D-03 



1.508857D+00 
2.5001 02C+00 
3.621 151D-03 
3. 762 1 2 ID-02 



EIG EKSY STEM CF OPTIMAL EEGCLATOR 



EULER-LAGRANGE EQUATIONS HAVE A REAL EIGENVALUE 
AT OR NEAE ZERO. 

C-LCCP OPTIMAL BEG. E- V ALU E S . . . DET (S I- F+G *C) . . 
0.0 : O.C ,- 1 . OOOOGD+OO: 

C-LCCP EIGHT EIGENVECTOR MATRIX M 



-7. CSS8C7D-C1 6.035 1 65D-0 1 
-7. 0823C7C-01 1.000CCCD+00 



C*i-1 



CONTROL EIGENVECTOR MAT BIX. 

- 1 .ill 176 1D+00 -1 . 5C4737D-02 

G ■*- *- *• ■£* ^ L C i X U — V* £ C TO R h xi a d X X • • U " 1 v • b 

-3. 592CS2D+CQ 2.167ES8D+00 
-2. 544E62E+G0 2.535E62D+00 

THE OPTIMAL FEEDBACK GAIN CONTROL MATRIX. . .C=BINV*GT*5 



5. lCS5i+GC -3.09671 + 00 

THE MCDAL CONTROL GAINS C*T. . 

-3 . C 5 66 9 6D+ 0 0 -5. 1C9451D +CG 

TEE CICSED LCOP DYNAMICS MATRIX F-G*C.. 

0.0 1.000CCCD+00 

4.1C9451D+C0 -3 .098696D+00 



ANALYSIS COMPLETE. CO YOU WANT ANOTHER RUN? 

TYPE "YES" 05 "NO". 

yes 

EC YOU WISH TO SiVZ THE "F"-MATRIX FROM THE LAST 
FUN TO EE USED IN THE FOLLOWING RUN? 

NOTE: THE M ATRI X WILL EE F EDIS PLAYED AT 
TEE PROPER INPUT SEQUENCE INTERVAL 
AND YOU WILL HAVE THE CPTI0N OF CHANGING 
INIIVIEUAL MATRIX ELEMENTS. 

TYPE "YES" OE "NO". 

yes 

DC YOU WISH TO SAVE THE "H"-MATRIX FROM THE LAST 
FUN TO EE USED IN THE FOLLOWING RUN? 

NOTE: THE MATRIX WILL EE REDISPLAYED AT 
TEE PROPER INPUT SEQUENCE INTERVAL 
AND YOU WILL HAVI THE OPTION OF CHANGING 
INDIVIDUAL MATRIX ELEMENTS. 

TYPE "YES" OR "NO". 

yes 

DC YOG WISH TO SAVE THE "G"-MATRIX FROM THE LAST 
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RUN 1C EE USED IK THE FCIICWING RUN? 

SCI E: IHZ MATRIX WILL EE FED IS FLAYED AT 
THE FRCPER INPUT SEQUENCE INTERVAL 
ARE YCU WILL HAVE THE OPTION OF CHANGING 
INDIVIDUAL MATRIX ELEMENTS. 

TYPE "YES" OB "NO". 



yss 



Author's rote: Since the s&me program options are to 

be run again, with only a change in or= 
of the cost matrix element values, the 
terminal output was deleted up to the 
point where program calculations resume 
in order tc avoid redundancy. 

ORDEB 01 SYSTEM = 2 

NUMEIE CF CONTROLS = 1 

NUMEE5 CF CESERVATICNS = 2 

NUMBER CF PROCESS NCISE SOURCES = 0 

THE SYSTEM MATRIX (»F"-MATRIX) ... 

C.O 1.00000 

-1.CC0CC 0.0 



EC YCU WISH TO CHANGE THE VALUE OF ANY MATRIX ELEMENT? 
TYPE "YES" CR "NO", 
no 

OPEN LC CP DYNAMICS MATRIX F. . 

0.0 0.100CC+01 

- 0. 1 C CCE + Q 1 0.0 



CPFN LC CP EIGENVALUES DET(SI-F).. 

0.0 , 1 . C000 CI+00: 

OPEN LC CP RIGHT EIGENVECTOR MATRIX T 

0. 0 -1 .000 CC0D + 00 

1. 0C00CCD+ CO 0.0 



64 



CEEF ICC? LEFT EIGENVECTOR MATRIX T-INV. . 

0- 0 1 .000 OCOD+OO 

-1. OOCCCCE+CG C.O 

TEE MEASUREMENT SCALING MATRIX ( "H" -MATRIX) 

0.0 0.0 

C..0 -1.0CCC0 

IC YCU WISH TO CHANGE THE VALUE OF ANY MATRIX ELEMENT? 
TYPE "YES" CP "NO", 
no 

MEASCSEMENT SCALING MATRIX H. . 

0.0 0.0 

O.C -0.100CE+01 

KCIAI MEASUREMENT SCALING MATRIX H (3 AH) *T. . 

C.C O.C 

- 1 . GCCOOOD+OO O.C 

ENIER THE OUTPUT MEASUREMENT COST MATRIX ( " A « -H ATEIX ) . 
II ME NS ION = # 0 ESERVAT ION S (NO) X # OBSERVATIONS (NO) 
THE ELEMENT A ( 1, 1) = 

C 

TEE ELEMENT A( 1, 2) = 



0 

THE ELEMENT A ( 2, 1) = 



C 

TEE ELEMENT A( 2, 2) = 

* 

4.000 1 

THE CUIPUT MEASCEEHENT CCST MATRIX (" A"-MATRIX) . . . 

C.O 0.0 

C.O 4.0CC10 

IC YCU WISH TO CHANGE THE VALUE 0? ANY MATRIX ELEMENT? 
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TYPE "YES' 1 C 3 "NO", 
r.e 

OUTPUT COST MATRIX A.. 

O.C 0.0 

0.0 Q.400CC+01 

THE CONTROL DISTRIBUTION MATRIX ( "G "-MATRIX) . . . 

0.0 

1 . ccooc 

rc YOU WISH TO CHANGE THE VALUE OF ANY MATRIX ELEMENT 
I YPE "YES" CB "NO", 
no 

ENTER THE CCN'TSCI COST WEIGHTING MATRIX ("5"-MATRIX) 
DIMENSION = * CON TROIS (NC) X # CONTROLS (NC) 

THE ELEMENT B ( 1, 1) = 

7 

1 

TEE CONTROL COST MATRIX B... 

1. CCOOO 

IC YCU WISH TO CHANGE THE VALUE OF ANY MATRIX EIEMEN1 
TYPE "YES" CB "NO", 
no 

THE CONTROL DISTRIBUTION MATRIX G. . 

0. C 

0. 1 CCOD + O 1 

KCIAI CONTROL DISTRIBUTION MATRIX TI*G. . 

1 . CCCOCCD+OO 

0.0 

THE CONTROL COST EATRIX B. . 

0. 1CC0E + 0 1 

CEEN LO CP TRANSFER FUNCTIONS... 

TF FCF INPUT NO. 1 AND OUTPUT NO. 1; 

NC FINITE ZEROS. Tr GAIN = 0.0 
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RESIDUES AT THE PCXES 



BEST DU 



5 



POLES 

REAL (A) IHAG(B) 

( O.C )+G( 1. COOOCC} ( 0.0 ) EXP (A*T) *C03 (E*I) 

( O.C ) *J ( -1.C000CC) ( 0.0 ) EX? (A*T) *SIN(E*T) 

TF FCF INPUT NO. 1 AND OUTPUT NO. 2: 



ORDER CF NUMER ATCR = 1 TF GAIN = -0.1000D + 01 



NUMERATOR EIGENVALUES (INCLUDING EXTRANEOUS ZERO VALUES): 
( O.C ) + J ( 0.0 ) 

( O.C ) +J ( 0.0 ) 



RZSIIUES AT THE POLES: 
E C L E S 

REAL (A) IN AG ( E) 

( O.C )+J { 1. COOOCC) 

( O.C )+J ( -1. COOOCC) 



RESIDUES 



( -1.000000) EXP(A*T) *COS (B*T) 

( 0.0 ) EX? (A*T) *SIN (B*T) 



EIGENST STEM CF OPTIMAL EEGULATOR 

C-LCCP OPTIMAL REG. E- V ALU ES . . . DET (S I- F + G *C) . . 

-1. 0C5C1E+CC:-9 . 950 1 PD-0 1: 

C-LCCE FIGHT EIGENVECTOR MATRIX M 

7.053368D-01 -7 .0 88 7 23D-0 1 
-7. 0 6 6723D- Cl 7 .053368D-01 

CCKTECL EIGENVECTOR MATRIX C*M. . 

1-41776 2D+ 0 0 -1 . 410691D+00 

C-LCCP CPT . REG- LEFT E-VECIOR MA TRI X. . M- IN V . . 
-1.41C6S1D+C2 -1 .4 17762D+02 
-1. 4 17762D+C2 - 1 .4 10 6 S 1D + 02 
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THE OPTIMAL FEEEEACK GAIN CCNTHCL MATRIX 
-3.1S74C-14 -2 .00 OCE + OO 



,C=3INV*G' 



TEE MODAL CONTROL GAINS - C*T. . 

-2.CCCC25D+00 3.1S7442D-14 

TEE CLOSED LOOP DYNAMICS MATRIX F-G*C.. 

0. 0 1 .0 00 C C OD + OO 

-1. 0CC00CC+C0 -2.000 C 25D + 00 



ANALYSIS COMPLETE. CO YOU WANT ANOTHER RUN? 
TYPE "YES" OR "NO". 



EC 



CPTSYSX IS NOW TERMINATED 

R; T = C.. 63/2 .60 23:33 :07 
rsccxc cff 

END RECORDING OF TERMINAL SESSION 



1 
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v. co^cinsicjss ahc be com men pat ions 
a. ccnciosioss 

Although originally developed fox the quadratic 
synthesis cf controllers for rotary-wing VTOL aircraft, the 
extensive modifications and enhancements of Hall's original 
work, coupled with its efficient and accurate eigensystsm 
solution routine, represent a powerful tool in the design of 
optimally controlled systems. 

In its present interactive form, OPTSYSX has been trans- 
ferred from the arena of high-level applied mathematics and 
numerical analysis tc the level cf control system engineers 
and students. It now represents an even more powerful 
educational tool, able to rapidly and effectively unlcck 
many n isunderstcod linear systems mathematical relation- 
ships . 

As an ultimate evaluation cf the computational abilities 
cf OETSYSX, the program was tested using an 82 X 82 matrix 
cf aircraft longitudinal motion equations for the VX-29 
experimental Fighter aircraft derivative, provided by 
K AS A - 1 ewar ds. 

For a system of equations of this magnitude, all program 
arrays were re- dime rsicned (as shown in Appendix A) , and a 
2-Megabyte virtual machine size was required. This system 
was run through the Medal Analysis option of OPTSYSX, 
requiring less than 90 seconds tc lead the system and 
complete all open-lccp and modal analysis calculations! 

Erogram results exhibited perfect eigenvalue correlation 
with these obtained from the John Edwards Control Ercgram. 
Additionally, OETSYSX provided complete longitudinal medal 
analysis, previously unavailable on a system cf this size. 
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use c 



It is hoped zhat 



£ th 



o 



interact 



V5 



* r* ■* • 
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version 
■will cx 

ccnxxcl 



fcill be encouraged; and that its expanded 
mulcts icxh interes tin and research on has 
problems, as *ell as more advanced designs 



a b 



T 

J. 

T» 



E. EiCCEMEHDATICHS 

Based cr. the resulzs of zh is tmesis, four arses amerced 
as possibilities for further research and szudy: 

1 . Proo ram Availab ilit y 

lbs use of CPTSYSX and similar design programs 
should be encouraged in all undergraduate and graduate level 
courses involved in the analysis and design of control 
systems. Toward this end, iz is recommended that CPTSYSX be 
placed ir the ncn-I2SL library of subroutines, makirg it 
easilv available to all potential users. 

2. Cc m cute r Gr ap hi cs 

The addition cf graphical plotting routines tc the 
program in the time and frequency domain would make CE1SYSX 
an even mere powerful tool in the design of many optimally 
controlled systems. 

3. F c rxhe r Modifications 

The present version cf the program should be modi- 
fied to include the CFTSYS 5 derivative input term improve- 
ments cf Liu [Eef. 3], and program sequencing during optimal 
filter synthesis should be examined. Various test runs 
indicate an area of conflict in that the program appears tc 
require the design of an optimal regulator prior to 
perfenirg any filter calculations. 
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APPENDIX A 

OPTSYSX PROGRAM LISTING 





C** ********************************4 ****** ***44******?********** ******** 

£** ********** ************ 
c********** OPTSYSX ********** 

C** ****** EY JOHN G. HCCEN ******** 

THIS PROGRAM IS A COMP IETELY INTERACTIVE ******** 

OETIMAL SYSTEMS CONTROL DSSIG S/STNTHESI3 ******** 

c **»,*,** PBCGR AH CAPABLE OF HANDLING 7 EH Y LARGE (30X80)+ ***»**•« 

C******** MCLTIVARIABLE SYSTEMS CF LIP EAR EQUATIONS. 

C** ****** 

c... ...*,., VERSION 1.8 11 SAH 1984 ********** 

i************ 
t ****** ************ 

IMPLICIT 8EAL*8 <A-H,C-Z) 

INTEGER IANS , ICL , IQ , 13 , IS S, IM , ITF 1 , 1 IF2 , 1 IF 3 ,1 FDFW , IE, IDESDG , IS ET, 

' 1 1PSD ,1 YtJ,INORM ,NS,NC,NOB , NG ,1 BEG , IDS T AE , IS ET , NROW , NCOL , IS A F, ISAH ,1 
2S AG, IG AM 

C LA5GE QBEEB SYSTEM (82 X 82) DIMENSIONS. 

C 

C 

c 

C *92 1 j&2 ,82) ^X ( 164, 164) ,GX <82 ,32)\HO (4 1, 62) ,31 ('l 64),D2]i64)/, 

C * EM (164 , 1 64) , Q (41 ,41) , GAM (82 ,4 1) , 3 NOR M (82,82) . 9 NOR3I (82 , 82) , 

C * DESTAB (8 2) , A A (8 2,82) ,BM (8 2, 4 1) ,CM (41 ,82) ,3 (4 1 ,4 1 ) ,3 STORE (32 .32) , 

C *JCF (16 4 ) RES (164) ,AY (82 ,3 2) , 3 c ( 164) , CC ( 164) , C? (32) , GW ( 16 4 , 4 1 ) , 

C * G V (164 ,41) , 3 Y (41, 164) ,HU (41 , 164) ,PBTT (16, 16) , DUM (82,85) 

C STANDABD FROGBAH DIMENSIONS. 

DIMENSION 
11(32) ^CV 

3 1 G N ( i 2 Jl) ,80 (32 ,32) .61 (6 4) ", 3 2 (64) , RM (6 4, 64) ,'Q (3 2,3 2) , GAil (32 ,32) ,9 

4 HCBM(3 2.32) , WNCBMI (32,32) ,3ESTAB (32) , A A (32 , 32) , BM (3 2 ,3 2) , CM (3 2 , 3 2) 
5,3(32,32) ,DSTCRE(32,32) , JCF (6 4) , RES ( 64 ) , AY (3 2, 32) , BB (64) , CC (6 4) , C? 
6 (32) , G 3 (64,64) ,GV (6 4,64) , HY 64,04), HU (64, 6 4) , PRTT ( 1 6, 1 6) , DUM (3 2,32 

H 

EQUIVALENCE ( 9 1 1 ( 1 , 1) ,G3 ( 1 ,~ ) ) , (Wll (1,1) , GV ( 1 , 1 ) ) , (9 2 1 ( 1 , 1 ) , HY ( 1 

i 

COMMON /EBOG/ IOL ,IQ,I3 f ISS , IF. ,ITF1 , IT? 2 , ITF3 , IPDF9 , IE , ID5T AE , I C EB 
1UG,ISST,IBEG,IPSC,IYO,INOBM 

DATA I Y / ' Y ' / , I Z / ' S'/ 

C SUPPRESS INDIVIDUAL 0NDERFLC9, 0VEEFL03 , DIVIDE CHECK, AND DECIMAL = 
C CCNVEET 2H8CB MESSAGES; PROVIDE SUMMARY OF ERRORS ONLY. = 

CALL ERRSET ( 2C 7 , 256 ,- 1 , 1 , 1 , 2 C 9) 

CALL ERRSET ( 2 1 5 ,256 ,-1 , 1 )_ 

C INITIALIZE PLAGS. 

I SAF=0 
IS AG=0 
ISAH=0 
IGAM=0 

TO CCNTINOS 
IBET=0 
ICL=0 
IQ=0 
18=0 
ISS=0 
IM=0 
ITF1 =0 
IIF2=0 
IT? 3=0 
I ?BF9= 0 
I E=0 

IDSTA3=0 
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0 

30 

40 

c 

5 0 

60 

70 

C— - 
30 

90 
100 
c — 



c 

110 

120 

130 

C 

c 

140 

150 

160 



I DEB UG -0 
ISET=0 
I F SD=0 
I 7 0=0 
INOR»=0 
IHEG=0 
N S=0 
NC=0 

ncb=o 

NG-=Q 

SCRN1 - 

CALL F ET CSS («CL2SCS3J » ) 

WRITS (5 .890) 

CALL 3 DC HAH (IANS) 

IF ( (I ANS.NE. IY) . AND. (IANS. NE.IZ) ) GO TC 30 
GC TO 40 
WHITS (5 , 88 0) 

GC TO 20 
CONTINUE 

I? (IANS. EQ. 12) GC TC 560 

SCHN2- 

CALL FRTC5S ( 1 CLRSC3 N •) 

WHITE (5 ,50 0) 

CALL RDCHAR (IANS) 

IF ((IANS. NE.IY) .AND. (IANS. NE.IZ) ) GC TO 60 
GO TO 70 
WRITS (5 ,880) 

GC TO 50 
CCNTIN UE 

IF (IANS. EQ. 12) GC TC 560 

1 SET-- 

CALL FSTC2S ( 1 CLRSCE N ') 

WHITE (5,910) 

CALL RDCHAR (IANS) 

IF ( (I ANS.NE.IY) . AND. (IANS. NE.IZ) ) GC TC 90 
GC TO IOC 
WHITE (5,880) 

GC TO 80 
CONTINUE 

IF (IANS.EQ.IY) I SET = 1 

CALL ? HTCSS ( * CLRSC E N 1 ) 

WHITE (5 ,570) 

CALL RDINT (IANS) 

ICL-IA NS- 1 

IF (10 L. EQ* 2) GO TO 350 

IQ 

CALL FBTCFS ( 1 CL3SCS N *) 

WHITS (5 .580) 

CALL RDCHAR (IANS) 

IF ( (IANS-NE.IY) .AND. (IANS. NE.IZ) ) GO TC 120 
GO TO 1 3 C 
WHITE (5.880) 

GO TO 110 
CONTINUE 

IF (IANS. EQ. 17) IC=1 
IF (IANS.EQ.IZ) IC=0 
IF (IOL.EQ.3) GO TO 200 

— — — — 1 R- 

CALL F3TCSS (’CLRSCSN ') 

W BITS (5,590) 

CALL RDINT (IANS) 

I3=IANS-J __ s 

CALL FHTCKS ( ' CLRSCB N ') 

WHITS (5,600) 

CALL RDCHAR (IANS) 

IF ( (IANS. NE. II) .AND. (IANS. NE.IZ) ) GC TO 150 
GC TO 160 
SHITS (5 , 

GO TO 
CONTINUE 

IF (IANS.EQ.IY) ISS=1 

IF (IANS.ZQ.IZ) ISS=G_ __ __ 

WRITS (5,610) 



C 

17 0 



73 



no no 



180 

190 

200 
c — 



CALL 3DCHAR (IANS) 
15 ( (I ANS.NE. 17) J 
GO TO 190 
WHITS (5.880) 

GO TO 17C 
CONTINUE 



,ND. (IANS. NE.12) ) GO TO 130 



IF (IANS.2Q.ITj 



IF (IANS. EQ. 12 
CONTIN US 
IF (IOL.EQ.3) IM=1 



IS=1 

IE-0 



CALL FHTCHS ('CLRSCEN ') 
WRITE (5,620) 

CALL 3DINT (IANS) 

ITF1=I ANS-1 

IF (IOL.EC.3) GO 10 240 



-ITF1- 



CALL F BTC MS ('CLESCSN » ) 
WRITE (5,630) 

CALL 3 DIN I (IANS) 
ITF2=IANS-1 

IF (IOL.EQ.3) GO TO 240 



-ITF2- 



C 

210 

220 

230 



CALL PSTCBS ( * CL BSCS N ') 
WRITS (5,640) 

CALL R DINT (IANS) 

ITF3=I ANS-1 



-ITF3 — 



220 



C 



CALL r ST CSS ( 1 CL3SCEN • ) 

WRITS (5,650) 

CALL RDCHA2 (IANS) 

IF ( (I ANS.NE. IT) . AND. (IANS. NE.IZ) ) GC TC 
GO TO 230 
WRITE (5,880) 

GC TO 21 G 
CONTINUE 

IF (IANS- EQ • I T ) IF DF W=1 
IF IANS. EQ. 12 IFDFWsQ 



-IFDFW- 



-IE- 



C 

240 

250 

260 



C 

270 



280 

290 

300 
C 



CALL PETCSS {'CLRSCEN ') 
WRITE (5 ,660) 

CALL 3DREAL ( A NS S ) 

IE = IDI NT ( ANS R ) 

IF (IOL-EQ.3) GO TC 300 






CALL FRTC2S ( ' CL3SCE N ') 

WHITE (5 ,670) 

CALL 3DCHAR (IANS) 

IF ( (IANS. NE. IT) .AND. (IANS. NS. 12) ) GC TC 
GC TO 260 
WRITE (5.880) 

GO TO 240 
CONTINUE 

IF (IANS . EQ. IY) IDST AS= 1 
IF (IANS. EQ. 12) I CST AB= 0 

WRITE (5 .680) 

CALL RDCHAR (IANS) 

IF ( (I ANS.NE. IT) .AND. (IANS. NE. 12) ) GO TO 
GO TO 290 
WRITE (5.880) 

GC TO 270 
CONTINUE 

IF (IANS . EQ. 17) ItSBCG=1 
IF (IANS. EQ. 12) I EE B UG= 0 
CONTINUE 



250 

• — IDE3UG- 
280 



CALL F ST CBS ('CL2SCSN ') 

^ — 0 v 

(IANS) 



-I PSD- 



WRITE J5.690V 
CALL RDINT 



IPSD=I ANS 
IF (IPSD.EQ.3) IP SD= C 
IF (IPSD.EQ. 0) GO TO 310 



CALL FRTCBS (’CLRSCEN 
WRITE (5 ,700) 



') 



-no- 



74 



n to 



10 

320 

330 
34 0 



c 

350 



CALL RDINT (IANS) 

IYU=IAN5-1 

IHORfl 

CALL FRTCMS ( ' CLRSCE N ') 

SR ITS (5 , 820) 

CALL R DR HAL ( A NS 5 ) 

IHOR21* IDINT ( ANSR) 

IF (IQ L. EQ. J ) GC 10 350 

CALL FRTCKS (•CLRSCEN •) 

SHITS (5,710) 

CALL RDCHAR (IANS) 

IF ( (I ANS.NE.I Y) . AND. (IANS. NE.IZ) ) GC TC 330 
GC TO 340 
WRITS (5.880) 

GC TO 32 U 
CONTINUE 

IF ( I A NS . EQ. IY) IBEG=1 
IF (I A NS * EQ. 12) I HEG =0 



CALL FRTCMS ('CLRSCEN ») 
WHITS (5,720) 

CALL RDRHAL (ANSR) 

NS = IDI NT (ANSS) 

IF (IOL-HQ.2) GO TO 360 



-NS- 



C — 



360 
C — - 



WHITE (5,730) 

CALL SDREAL (ANSR) 
NC-IDI NT (ANSR) 

WRITE (5 ,740) 

CALL RDRHAL (ANSR) 
NC3= ID INT (ANSR) 









-NG- 



W HIT S (5,750) 

CALL RDREAL (ANSR) 
NG=IDINT (ANSR) 
CCNTINUE 



•SLAG SETTINGS- 



WRITS 


(6 ,760) 




WRITE 


6 ,770 




WRITE 


6 ,780 


IOL , 


1 , ILST A B 




WRITE 


(6 ,790) 




WHITE 


6 ,800 


IPS c 


WRITS 


6 ,810 


ES, hi 



C 

370 



380 

390 



BEGIN CALCULATIONS 

N2=2*NS 

CALL INN HE ( NS , NC ,NO E ,N G , N2 , ACI, B , BA ,CI ,CR ,CQ, CWI ,C WR , D, FBGC , FBG H , 
1G,GAM,GM,GN,HC,D1 ,D2 .PRO. EM ,2C ,Q,3C, HR, WI ,W1 1 , W2 1 , X , WNCRM, S NO FBI ,D 
2ESTAB, A A , 3M , CM ,JCF,RHS, AY ,33 , CC ,C? , G S ,G7 , fl Y , HU , DSTO R E, IS AS , IS AH , IS 
3 AG , IGA M , IHET, ERTT,NRCW, NCCL) 

1 RET 



C 

400 



410 

4 20 
430 



WRITS (5.830) 

CALL RDCHAR (IANS) 

IF ( (IANS. NE.IY) .AND. (IANS. NE.IZ) ) GC TC 380 
GC TO 390 
WRITE (5,880) 

GC TO 370 
CCNTINUE 

IF (IANS.EQ.IY) GC TC 400 
IF (IANS. EQ. 12) GC TC 560 

2 i ISAF- 

CQNTINUE 
IF (IRET.EQ. 1) 

I? (ISET.EQ.1) 

CALL FSTCMS ( * ( 

WRITS (5 .840) 

CALL RDCHAR (IANS) 

IF ((IANS. NE.IY) .AND. (IANS. NE.IZ) ) GO TC 4 20 
GC TO 430 
WHITE (5 ,880) 

GC TO 4 1 C 
CCNTINUE 

IF (IA NS. EQ. IY) I SA F “ 1 



GO TO 10 
GO TO 10 
' CLRSCEN •) 



75 



c 

440 

450 

460 



470 

C— 



IF (IA 



480 

490 

500 



51 0 
C— 



520 

530 

540 

550 



C 

560 

C 

570 



IF (NO 
CALL F 

'« p T Tt? 

CALLER 
IF ((I 
GC TO 
KEITH 
GC TO 
CCNTIN 
IF (IA 
IF (IA 
CCNTIN 



NS.EQ.IZ) I S AF = Q 

ISAH-- 

3-EC.O) GC TO 470 
RTCiiS ( * CL2SC5 N •) 

(5,850)' 

DCnAR (IANS) 

ANS.NS, I Y) . AND. (IANS. N2.IZ) ) GC TC 450 
4 o C 

(5 .580) 

44 C 
OE 



NS.HQ.IY) IS AH = 1 
NS.EQ.f' ' 



OE 



IZ) ISAH=0 



IF (NC 
CALL F 
KEITH 
CALL E 
IP (d 
GC TO 
KEITH 
GO TO 
CONTIS 
I? (IA 
IF (IA 
CONTIN 



IF (NG 
CALL ? 
KEITH 
CALL R 
IF ((I 
GO TO 
KEITH 
GC TO 
CCNTIN 
IF (IA 
IF ]lA 
CCNTIN 
GC TO 



ISAG- 

* E C - 0) GC TC 510 
RTC3S ( • CLSSCBN * ) 

(5 ,860) 

DC EAR (IANS) 

AMS.NE. IY) . AND. (IANS.ME.IZ) ) GC TC 490 
50C 

(5 .880) 

480 

UE 

NS.2Q.IY) I SAG= 1 
NS.EQ.IZ) IS A G=0 
UE 

IGAM- 

.EC. 3) GC TC 550 
ET C!*.S ( • CLRSCEN • ) 



(5 .870) 

DCfcAE (IANS) 

ANS.NE.IY) .AND. (IANS. NE.IZ) ) GC TC 530 
540 

^S e80 > 

OE 

NS.EQ.IY) IG A M = 1 
NS.EQ.IZ) IGAM=0 
OE 
10 



W BITE 
STOP 



(5,920) 



-T E3MIN ATE 



FORMAT f25X,24HGEcSERAL OPTSYSX CPTIO NS : ,// , 1 OX ,3 5 HO PTION 1 — SYST 
12S ANALYSIS 4 ITHOOT , / , 2 2 X , 3 5 HCEEN-LO CP EIG £n S Y STEfl CALCOL AT IONS . ,/ 
2/. 10X, 42HCPTICN 2 “ SYSTEM ANALYSIS WITH OP EN -LOO? , /, 22X, 2 5 HEIGEN 
3SYST2M CALCULATIONS. ,// , 10X,39HCETION 3 — OPEN-LOOr 2IG2NSYSTEM F 



o . d nciuui.cij.uci iuucci iCr/ti. 

7 , 25HOR STEADY-STATE ANALYSIS. ,// , 1 5X , 30HSELECT Ail OPTION: 1,2,3, 0 

8R 4.) 

580 FORMAT (//,5X,46HCO YOO DESIRE RMS VALUES OF STATE AND CONTROL?,// 
1 , 10X.1 9HTIPS "YES" CR "NO".) 

590 FORMAT (/ , 20X , 30 HCFTSYS X LQ R/CLASSIC AL OPTIO NS : , //, 1 0X ,43 HOPTION 1 
1 — OPTIMAL FILTER AND/OR REG OLATCR, / , 22 X, 37 HS YMTHES IS WITH NC EXT 
2EBNAL "C" OR "K",/,22X, 13HMATSIX INP OT . , // , 1 OX ,4 3 HO PTION 2 — OPTI 
3MAI FILTER AND/OR RE GOL AT CR, / , 22 X ,27 HS Y NT H ES IS WITH EXTERNAL "C",/ 
4 ,22X, 1 3HMATRIX IN EOT . ,// , 10 X , 5 3HCPTI CN 3 — OPTIMAL FILTER AND/OR 
5REGULATOR,/,22X,27HSYNTHESIS WITH EXTERNAL "K" ,/ , 22X , 1 3HM ATR IX INP 
6 0 T . , // , 10x,43fiCllICN 4 -- OPTIMAL FILTER AND/OR REGU LATOR , / , 22X , 35 
7 8S YHTH ESIS WITH EXTERNAL "C" AND "K" ,/, 22X, 1 3H MATRIX INPUT. ,//, 1 CX 
8, 32 HS ELECT AN OPTION: 1, 2, 3, OR 4.) 

600 FORMAT (//,5X,SCHEO YOO 'WISH TC DETERMINE THE STEADY-STATE RESPONS 
1E,/,3X,27EFOR A CONSTANT CISTURB ANC E ?, // , 1 OX , 1 9HT YPE "YES" CR "NC" 
2.) 

610 FORMAT (5X,U7HEO YOU WISH TO EETERMINE THE MODAL DISTR IBUTICN ,/ , 8X 
1 , 1 8H AN D GAIN 2 ATR ICES?,// , 1 0 X , 19HTY P E "YES" OR "NO".) 

620 FORMAT (//, 5X , 36 HCPE N-LO OP TRANSFER FUNCTION OPTIONS :,//, 1 OX , 53HCP 
1TION 1 — NO 0 EE N-LO CP TRANSFER FUNCTIONS COMP UTED. ,// , 1 0 X , 48 BC ETI 
20 N 2 — POLES, RESIDUES , AND ZEROS COMPUTED. ,//, 1 0X , 42 HOPT ICN 3 — 
3 ONLY POLES ANC ZEROS COMPUTED. ,//, 1 OX, 45 HCPTION 4 — ONLY POLES A 
UND RES ID EES CC MPCTED . ,// , 10 X , 3 2HS ELECT AN OPTION: 1, 2, 3, OR 4.) 

630 FORMAT ( //, 5 X , 32H NCI SE TRANSFER FUNCTION OPTIONS : , // , 1 OX , 4 9 HOPTICN 
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640 

650 

660 

670 

580 

690 

700 

710 

720 

730 

740 

750 

760 

77 0 

780 

790 

800 

810 

820 



1 1 -- SO NOISE TRANSFER FONCTICNS CO MPUTSD. r // , 1 0X, 48HOPTIOU 2 — 

2ICIES , RESIDUES, AN C ZEROS CONFUTED. ,// , 1 OX, 42HQPTION 3 — CHIT FO 
3LES AND ZEROS CONFUTED. ,//, 1 0 X ,45 HOPTIQ N 4 — ON LY POLES All L RESID 
4UES COMPUTED, , // , 10 X , 32 HS EL ECT AN OPTION: 1, 2, 3, OR 4.) 

FORMAT f // , 5 X , 3 8HCC M PENS A TO R TRANSFER FUNCTION OPTIONS :,//, 1 0 X , 4 9H 
1CFTION 1 — SC CCMP. TRANSFER FtJNCTI C NS COMPUT ED .,//, 1 OX , 4 8 HOFTTCN 

2 2 — POLES, RESICUES, AND ZEROS COMPUTED. ,//, 1 0 X , 4 2 HO PTIO N 3 — O 

3NLT POLES AND ZEROS COMPUTED. ,//, 10X ,45HCPTION 4 — ONLX POLES AND 
4 RESIDUES COMPUTED. ,//, 15X, 45H NOTE: A COMPENSATOR TRANSFER FUNCTI 

50 N CAN 3E ,/,22X. 33HCCMPUTZD O MY IF EOT B A REGULATOR, / . 22 X , 26 HA N E 
6FILTZR ARE SlNTE IS1ZED, / r 22X , 14HAND/CR INPUT. ,//, 1 0 X , 32 HS ELECT AN 
7 OPTION: 1, 2 , 3, C3 4.) 

u TriMfn / / / c v IQ'JUTTT » E'C'-rs ronn » on n t rn n t o rrm -m *t u i-norv / c •*> C f i 



1 



FCSMAT (//,5 X , 39 HWIL I A FEED-FCRWARD DISTRIBUTION a ATR IX , / , 5X , 25 H 
1"E" - MATRIX RE INPUT ? , // , 1 5X , 1 9HT Y PE "YES" 03 "NO".) 

FORMAT (/,5X,63H THIS OPTION DETERMINES THE CRITERIA FOR DECIDING 

WHEN A MARKOV , / , d X , 5 3H? A R AM STE a IS ZERO-THE 3ARKOV PARAMETER INDIC 
2ATES THE CRD E R , / . 8X , 54H O F THE NUMERATOR POLY NOMIAL OF EACH TRANSFE 
3 R FUNC TICS. ,//,3X,52f.ALL "N" 2 EROS OF THIS POL YNOHI A L ARE PRINTED 
400T AND, / .3X, 52HTHIS TEST TSLIS HCN 3 ANY EXTRA ROOTS EXIST AT Z = 
50. ,/,3 X, 41HL2SS THAN 10.0** -IE IS CONSIDERED ZERO . ,//, 3 X , 4 7H THE 
6 DEFAULT VALUE CF THIS PARAMETER IE IS 6. ,/, dX , 28HIN OTHER '-.'ORES 
7, IS = 1 .CE-o. ,//,10X,66HI? YCG DESIEE A DIFFERENT BARKOV CRITERIA 
9, TYPE TEE INTEGER V ALU S. ,// , 1C2 , 48H IF YCU DESIRE THE DEFAULT VALU 
9 E , TYPE "0” ZERO ) 

FORMAT (//,5X,61HCC YOU DESIRE TO SYNTHESIZE A STABLE FILTER OR R 
1 EGULATOR 3 Y , / ,3 X ,3« HDEST AB ILIZING THE ORIGINAL SYSTE3?,// , 1 2X, 52H 
2 NOTE: WCSKS FCR FILTER OR REGULATOR EUT NOT FOR BOD H,/ ,20 X , 1 7 HIN T 
3HE SAHE SON. ,// , 10 X , 19 HT YP E "YES" OR "NO".) 

F0R3AT (52,53 HDO YCU DESIRE TC PRINT THE EULER -LAGR ANG E EIGENSY3TE 
' 1 3,/,3X ,50H?RICE TC D ECO 3POS IT ICN FOR CHECKING THE PROGRA3 ?,//,10 
2 X , 15HT YP E "YES" CS "NC".) 

PC 33 AT (//,5X , 39HFOWER SPECTRAL DENSITY PSD OPTION 1 :,//,10X,53 
1HCPTION 1 — CC3PCTE THE PSD C? THE OUTPUTS AND/OR r HE ,/ , 2 2 X , 4 9 H CO 
2 NT SOLS OF THE CONTROLLED SYSTEM WHEN FORCED 3Y ,/ , 22 X ,4 5HPP.OCE5S AN 
3D SSASCREHEST NOISE. NOTE: ECTH A, /, 2 2 X , 46 HR EGULATOR AND A FILTE 

UR 3UST BE RESIDENT IS TH E ,/ , 2 IX, 23HP BOG S A3 TO USE THIS OPTION. , // 
5 , 1 CX, 5 3 HC ETIO N 2 -- SABE AS OETICN 1 ABOVE 3UT ONLY PRINT THE,/, 22 
6X,34HRESIDUES CF EACH TRANSFER FUNCT ION ,/ ,22X, 28HUS ED IN THE PSD C 
70BPDTATI0N. ,// , 1 0 X, 2 4H0PTI0N 3 — NOT DESIRED. ,///, 1 0X ,29 HSEL ECT A 
8 N OPTION : 1, 2. OS 3.) 

FORMAT ( / / , 5 X , 39 H POWER SPECTRAL DENSITY PSD OPTION 2 :,//,10X,35 
1HCPTION 1 — PSD CUT PUT NCT DESIRED. ,//, 1 0 X , 33 HOPTION 2 — COMPUTE 

2 ONLY OUTPUT P S C° S. , //, 1 0 X, 3 9 6CPTION 3 -- COMPUTE ONLY CONTROL PSD 

3 0 S . - // , 1 CX,50HCPTICN 4 — COMEUTE BOTH OUTPUT AND CONTROL PSD°S.,/ 
4/, 15X, 32ESELECT AN OPTION: 1, 2, 3, OR 4.) 

FORMAT ( // ,5 X , 39H DC YCU DESIRE REGULATOR SYNTHESIS ONL Y? , // , 10X , 19 
1HTYPE "YES" OR "NC".,//) 

•47BENTER THE * OF STATES NS OP THE SYSTEM MATRIX,/, 
CF CONTROLS NC OF THE CONTROL SYSTEM 

OF PROCESS NCISE SOURCES NG OF THE,/ 

DTE R SETTINGS FCR THIS RUN ABE AS FOLLOWS: 

,2HIM,2X,4HITF1,2X,4HI 
,2X,4HISET,2X,6HIESTA8 



FORMAT f/,5X,47EENT: 

1 5X , 1 3H "F"- MATRIX .) 

FORMAT (/.5X.56BENTER THE # 

1 MODEL / 5x 1 3 H "G"— MATRIX 

FORMAT'’ (/ ,52, 54BESTEH THE #‘oF MEASUREMENTS OR OBSERVATIONS NO O 
1 FTHE ,/ , 5 X , 1 3 H "H"-MATRIX .) 

FORMAT ( / ,5X , 4 8HENTE R THE # 

1 , 5X , 1 7 H "GAMMA "-MATRIX .) 

FORMAT (5X,52H£LAG/FABAME1 
1 ,/) 

FORMAT (1X,3HIOI,2X, 2HI Q , 2X,2HIR,2X, 3HISS,2X,; 

1TF2, 2X ,4HITF3 ,2X,5HIEDFW,2X,2EIE,2X,6HICEBUG,; 

Format (U,i2,3x,i2,3x, I2,2x,i2,3x,i2,3i,i2,ux,i2,4x,i2,4x,i2,4x,i 
1 2 3 X X 2 61 12 5x X2 /) 

f6r 3AT ( 1X*4 Hfis6,2f ,3HIIU, 2 X , 5HINOR 3 , 2X , 4 HIR EG , 21 , 2 HNS , 2X , 2H NC , 2X 
1 , 3HN03 ,2 X ,2HNG ,/) 

FORMAT (2X,I2,3X,I2,4X,I2,5X,I2,3X,I2,2X,I2,3X,I2,2X,I2,//) 

FORMAT (22, ~~ 

1 3 ,// , 2 X 



( 2 2 , 1 7 HCRDES CF SYSTEM = , 13 , // , 2X , 20H NU M 8 ER OF CONTROLS = ,1 
_ ,24HNUMEER OF OBSERVATICNS = , I 3 , //, 2X , 3 3H NU MB ER OF PROCESS 
2NCISE SOURCES =,I3,////) 

FORMAT (SX,53HDETERHINE THE NC EM ALIZ ATICN PARAMETER INORM FCR TH 
1E,/,5X,d5HPCWEE SPECTRAL DENSITY PSD OPTION YOU HAVE PREVIOUSLY. 
2/,3X,52HCHOSEN. TVC PSD NCRMA LIZ ATIO N METHODS ARE A V AILA8 L 2: , // . 10 
3 X , 54 HM ET HCD 1 — FSD IS NCRMAIIZED BY THE I-NORM°TH PROCESS ./, 2 1 X , 
429HNOISE MINUS "Q" INORM, INORM .,/,21X,49H NOTE: "Q" IS AN OPTIMAL 

5 STATE WEIGHTING MATRIX. ,/,21X,34HIN THIS METHOD, INORM =0,1,..., 
6NG.,// ,1 CX,63HMETKOD 2 -- PSD IS NORMALIZED BY THE INORM -NG °TH 
7MEASUREMENT,/,21X,39HNOISE MINUS "R" INCRM - NG, INORM - NG . ,/,21X. 
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3,51ri NOTE: "H" IS AN OPTIMAL CONTROL WEIGHTING MATRIX. ,/,2lX,44HI 
9N THIS METHOD , I NCR M = NG + 1,...,NG * NOB., // f 10X,51HSZLECT AN IN 

$T£GEH FROM 0-16 REPRESENTING YOUR ESD ./, 15X, 27HNO REALIZATION R £Q 
S UIREME NTS . # // # 1 OX ,5 3 HI F PSD NORMALIZATION 13 NOT DESIRED ENTER "C" 
3 ZERO ) 

830 FORMAT j;5X,43HANAI.YSIS COMPLETE. DO YOD SANT ANOTHEB R UN? , / , 1 5X , 19 
1HTYPE "YES" OB "NC".) 

840 FORMAT (///, 5 X ,4 3EDC YOU WISH TO 3AV E THE """-MATRIX PROM THE LAST 

1 ,/,5X, 36 £sUS TC BE USED IN THE FOLLOWING RUN ? , // , 5 X , 39 HNOT E : THE M 

2 A TBI X SILL 3E BEDIS PLAYED AT , / , 5 X , 34 HTH S PKOPEB INPUT SEQUENCE I NT 
3ERVAL,/,SX,40HANC YOU SILL HAVE THE OPTION OF CHANG ING ,/, 5 X , 2731 SD 
41 VIDUA L MATRIX E I E ME MS . , // , 1 5 X , 1 9HT Y P E "YES" 03 "NO".) 

850 FORMAT (///, 5X ,48EEC YOU SISH TC SAVE THE "H"- MATS IX FROM THE LAST 
1 , /,5X, 36HRUN IC £ E USED IN THE FOLLOWING BON? , //, 5 X ,39 HNOT E: THE M 
2 ATBIX SILL BE REDISPLAYED A? . / , 5X , 34 HTH E PROPER INPUT SEQUENCE I NT 
3ERVAL,/,5X,40EANC YOU SILL HAVE THE CFTION OF CH ANG ING ,/ , 5 X , 27HI ND 
41 VIDUA L MATRIX ELEMENTS. ,//, 15X.19HTYEE "YES" OR "NO".) 

860 FORMAT (/// , 5 X , 4 8 EDO YOU SISH TO SAVE THE "G"- MATRIX FROM THE LAST 
1 ,/,5x,3oHRUN TC EE USED IN THE FOLLOWING R UN ? , // , 5 X , 39 HNO T E: THE a 

2 AIRIX HILL BE REDISPLAYED AT . /, 5X , 34 HTH E PROPER INPUT SEQUENCE INT 
3EBVAL, /, 5X,40HAND YOU SILL HAVE THE OPTION OF CH ANG ING ,/ , 5 X , 27 HI ND 
41 VIDUA L MATRIX EXE ME STS . , // , 1 5 X , 19KT YPE "YES" OR "NO".) 

870 FORMAT (///, 5 X ,5 2ECO YOU SISH TO SAVE THE "GAMMA "-MATRIX FROM THE 
1LAST,/,5X,36HBUN TO EE USED IN THE FOLLOWING R UN ? , // , 5 X, 3 9 H NOTE : T 
2HE MATRIX HILL EE REDISPLAYED AT,/, 5 X , 3 UHTHE PROPER INPUT SEQUENCE 

3 INTERVAL, /,5X ,4CEANC TOO SILL HAVE THE CFTION OF CH AN GING ,/ , 5 X , 27 

4 HINDI'/ ID 0 AL MATRIX E IEMS NTS . , // , 1 5X, 19 HTY P E "YES" OR "NO".) 

880 FORMAT ( IX, 5 1 H WARNING: I.MPROPE5 DATA ENTRY! ENTER "YES" OR "NO".) 
890 FORMAT ( 5 X ,5 9 HOE I SYS X IS A COMPLETELY INTERACTIVE OPTIMAL SYSTEMS 
1 CONTROL,/,aX,55EPBOGFAM. IT Will SOLVE NUMEROUS CONTROL PROBLEMS 0 
2 N THE,/,8X,45 hECLICWISG TYPES CF SYSTEMS CONTROL SQU AT IONS : , // , 1 5X 
3 , 35HXDOT = ? *X ♦ G *0 + GAM * (S 4 WO) ,//, 20X , 22 HM EAS UREM ENT EQUA 
4T ION — ,//.15X,21H7 = H * X *■ C *W * V , // , 20 X , 29 HREGUL ATO S PERFORM 

5 A 11 C E INDtS — ,//,15X,42HJ = 1/2 * INTEGRAL (Y * A «Y + U * B *U)DT, 
6//.20X, 32 ESTATE FEEDBACK GAIN DEFINITION — ,//,25X,10HU = - C *X,// 
7,15X,45HDC YOU SIEH TO CONTINUE? TYPE "YES" OR "NO".) 

900 FORMAT (25X,14B — DATA ENTRY-- ,//, 5X , 49H ALTHOUGH OPTSYSX IS 3PECIFI 
1CALLY DESIGNED TC RE AD , / , 5X , 4 88 ALL MATRIX DATA INTERACTIVELY, SEVE 
2 S AL ALTERNATE, /,5X, 3 1HMETHODS ARE AVAILABLE TO USERS : , // , 1 0 X , 43 H ME 
3THCD 1 — THE "F","G" AND "GAMMA" MAT RIC ES , / , 1 3 X , 3 7 HM AY 3E READ FRO 
4 M SEPARATE DATA FILES. , // , 1 0 X , 5 OHMET HOD 2 — THE r 'F","G", AND "GAMMA 
5" MATRICES MAY BE ,/ , 13X , 4 5HEX £ IICITL Y DEFINED WITHIN SUBROUTINE "S 
6ETUP". ,//,10X,52H NOTE: IN EITEER CASE, THE USSR SHOULD OBTAIN A C 

70PY,/, 17X.34HCF TEE PROGRAM LISTING AND EX AM IN E, / , 1 7 X , 39HTHS EXAMP 

8 LES CONTAINED IN S/R "SETUP". ,// , 10 X , 4 5HDO YOU WISH TO CONTINUE? 

9 TYPE "YES" OB "NC".) 

910 FORMAT (//,5X,U6HDO YOU SISH TC INPUT TEE "G", AND "GAMMA",/, 

1 1CX,40HMATRICES FROM SUBROUTINE "SET OF" IAS TH E, / , 1 0 X , 40 HM ETHC D DS 
2SC5I3ED ON THE PREVIOUS SCH ES N ? ,//, 1 SX , 19HTY PE "YES" OR "NO".) 

920 FORMAT (//,41H OPTSYSX IS NOS TERMINATED ,//) 

END 
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non 



_5^52H TIN£ SET0F (£A,G,GAM, NS ,NC,NG) 

“HPLICIT'iilLiSTA^aTo^zr - 

DIMENSION 3A (NS, MS) ,G(NS,NC) ,GAM(NS, NG) , DU M ( 32 # 35) 

CCMMC N /f HOG/ IOL,IC ,IR , ISS ,1 t,ITF1 , ITF 5 , ITF 3 , IF DFW ,IE, IDSTAB,IDS3 
1 UG,ISE7, 1 BEG, IrSD , IY U ,INO RM 



FILE C£F I NITIONS 



CALL F BTC MS ( * FIL ED E F • * 03 1 . ' DISK »,'X29A82 

1 'DATA *, 'A *) 



C THIS IS AN EXAMPLE CF AS 82 X 95 CATA FILE X29A82 DATA A1 READ FROM 
C A USER'S DISK AND CONVERTED (FROM A "DUMMY'* ARRAY NAMED * DUM 1 ) TC A 
C SYMMETRIC ARRAY. THE FORMAT STATEMENT MUST MATCH YOUR DISK DATA 
C FORMAT OR THE PROGRAM JILL FAIL! NOTE: ALL PROGRAM DIMENSIONS 
C MUST EE ENLARGED ACCORDINGLY FOR A SYSTEM OF THIS SIZE. 

C 

READ (3,50) ( (LUM (I, J) , J= 1, 85) ,I = 1,N S) 

DO 20 1-1 , NS 
EC 10 J=1,NS 
BA (I, J) =DUM (I , J) 

10 CONTINUE 

20 CONTINUE 

C 

C THESE ARE EXAMPLES CF SEVERAL POSSIBLE METHODS OF ARRAY GENERATION 
C WITHIN SUERCUTINE SETUP. THE "GAM" ARRAY WAS SET TO ZERO SINCE NO 
C "NOISE” WAS PRESENT, AND THE NCN-ZZRC ELEMENTS OF THE "G" ARRAY WERE 
C EXPLICITLY DEFINED. THEY COULD ALSO BE READ FROM FILES AS ABOVE. 



C~ 

50 



DC 40 1=1, NS 
DO 30 J*1.SC 
GAM (I, J) =0. OD 00 



_ (I. J) =0 . CD-*- 0 0 

g Ja5, i) =c. ioooe+o i 

CONTIN UE 
CONTINUE 
RETURN 



FORMAT (5(212.4)) 
END 
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c 

c=- 



SUBROUTINE CHECK ( E P S , N C , NG , N C , IR ET) 

CONSlSTEiiCY OF REQUESTED OPTIONS. 



CKECXS 



C- 

c- 



10 
c- - 



D0U3LE PRECISION EPS 

COMMON /PROG/ ICl,IQ,IR, ISS ,IM,ITF1, ITF 2 , ITF 3 , IFD F W , 12 ,IDS T A B , I 0 EB 
1UG.ISST, IRSG f IPSD .IY(J,I NORM 

SET MODAL ANALYSIS 5 HEN OL EIGENSY S OS OL TF REQUESTED 

IF (I il .EC. 1 .ANC. IOL . EQ. 0) IOL= 1 
IF (ICL . EQ - 3 .OR. ITF 1 .NS. C) IM= 1 

CHECK TO SEE IF H MATRIX INPUT 

IF {NO • BE. 0 .Ofi. ICL . GE. 2) GO TO 10 
WRITE {5 , 90) 

IRZT=1 

RETURN 

CONTINUE 

TRANSFER FUNCTION CHECKS 

IF (IE . EC- 0) I£=6 

SPS = 10.** (-IS ) 

OPEN LOOP TF 

IF (IT F 1 .EQ. 0 .CR. NC . NE . 0) GO TO 20 

WRITE (5,100) 



C — 
20 



30 
C — 



c— 

40 



50 

c — 



60 

70 

80 



IHET=1 

RETURN 



IF (IT F3 .EQ. 
IF (IS EG .EQ . 
WRITE (5,110) 
IR ET= 1 
RETURN 
CONTINUE 



TF- 



CCM FENS ATO R 

0) GO TO 30 
C .AND. (NC . NE. 0 



.AND. NG . NE. J) ) GO TO 30 



NOISE 

0) GO ro 
.AND. NC . 



TF 

40 

NE. 0) 



IF (IT F2 .EQ. 

IF (NG . NE. 0 
W BITE (5,120) 

I HZT = 1 
RETURN 

DESTABILIZATION 

IF (IDSTAE .EQ. 0) GC TO 50 
IF NC .EQ. 0) GC TO 50 
IF (NG . NE. 0) IR EG= 1 



GO TO 40 



RESTRICTIONS* 



C 

90 

100 

110 

120 

130 

140 



WRITE (5,130) 

IF (IR EG .EQ. 1) 
IRET= 1 
RETURN 
CCNTIN UE 

PSD 

IF (I?SD .EQ. C) 
IF (IPS D .LT. 0 
IF (IYU .IT. 0 
IF IlNCRM .LT. 

GO TO 70 
WRITE (5,140) 
IRET-1 
RETURN 

I? (13 EG .EQ. 
WRITS (5,150) 
IRET- 1 
RETURN 
CONTINUE 
RETURN 



GO TO 50 



INPUT- 
GO TO 80 
OR. IP S D 



.OR. IYU . GT. 
0 .OR. IN 0 Rrt 



.GT. 3) GO TO 6 0 
3) GO TO 60 
. GT. NG + NO) GO TO 60 



0 .AND. NC .NE. 0) GO TO 30 



FORMAT (// ,5X , 
1 ,//) 

FORMAT {//, 5X , 
1 1 OX ,26HTC COMP 
F C 3 M A I (//, 5 X , 
1 / , 5X, 44H IN THE 
21 REG MUST = 0. 

FORMAT ( // , 5 X , 
1 £D ,/ , 5 X , 47HI. E 
FORMAT { //, 5 X , 
1 , 5X,33HOR FILT 
2. THE REGCLATC 
FORMAT { // , 5 X , 
1 ,//) 



4 9 E H - MATRIX MUST 3E INPUT, 1.3. "NO" MUST BE > 0. 

INPUT, 1.2. NC MUST BE > 0.,/, 



46 H (G) MATRIX MOST BE 
UTE OPEN LO CP T. F.,// 
48 BBEGULATO R ANC FILIE 
SAME RUN TC COMPUTE C 
; " NC" AND "NG" MUST 3 
51HNOISS T. F. CALCULA 
. I REG MUST = 1 . ; "NC" 
47HE2STABILIZATICN OPT 
ER HUT NOT BOTH SIMULT 
R OPTION IS AUTCMATICA 
49H ********** INCCNSI 



k SYNTHESIS MOST BE REQUESTED, 
CM PENS AT OR T. F . , / , 5 X , 4 7 HI . E . 

E > 0. ,//) 

TED ONLY WHEN REGULATOR DESIGN 
AND "NG" MUST 3E > 0. , //) 

ION DESIGNED FOR A REGULATOR,/ 
ANZOUSLY. ,//, 5X,55HIF "NG" > 0 
LLY SET!,//) 

STENT PSD INPUT FLAGS ******** 



80 



150 rOBNAT ( // , 5 X , 44 fiEOT E A REGULATOR AND FILTER «UST 3E RESIDENT,/, 10 
1X.42HTO CCiP'JTE TFE PSD Of A CONTROL LED SX STE3 ! , / , 1 0 X, 42 HI . 2. IREG 
2 MUST 3E 0. ANE "NC" .1UST BE > 0.,//) 

ENE 



I 
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SUBROUTINE INNER (NS , NC , NO, NG , N2 , AC L , B , 3 A , Cl , CR , CQ , CW I , CW 3 , D , F3GC , 
1 P3GZ ,G ,G AM,GM ,GN ,HO ,D1 , D2 ,FRO ,SM, 3C, Q,SC, WR, WI, W 1 1 , W21 ,X, NNOP M, WNO 
2BBI , DEST AE, AA, EM, CM , JCP r RE5,AY,33,CC ,CP ,G W ,G V, HI , HU , DS TORE , IS A? , 13 
BAH^ISAG^IGAM^ I HSI rPHT?^ NR CW , NCC1)_ 

IMPLICIT BEAL^A-H^i” 

C 




c- 

c- 

c- 

c- 

c- 

c- 

c- 

c- 

c- 

c- 



6 ,GV (N2 , II C) , BY (MO, II 2 ) , n r J (NC,N2) , DSTOB 2 (MS ,NS) ,?RTT (16,16) 



COMMON /PROG / ICL,IQ ,I2,ISS , IK, IT FI, IT F 2 , ITF3 , I? DFW , IE , IDS TAB r ID £3 
10G,ISET, I5EG, IPS C , IY U, IN 0 RM 

RZAL*4 F ST (20) 

OUTPUT OPTIONS 

■ ICL= 1 IF TEE OPEN ICC P EIG BUSY STEM IS DES 13 ED- -OTHERWISE IOL-O 

•10 = 1 IF THE RMS VALUES OF THE CC NT ROL - AN C STATE ARE TO 3E FOUND 

■IR=0 IF OPTIMAL FILTER AND REGULATOR SIGZNSYSTZMS ARE TO BE FOUND 

•IR=1 IF EXTERNAL C MATRIX IS SUPPLIED 

•IR-2 IF EXTERNAL K IS SUPPLIED 

* IR = 3 IF EXTERNAL C AND X A HE SUPPLIED 

■ISS=1 IF STEADY STATE VALUES ARE TO 3E DETERMINED 

■ IM= 1 IF MODAL STATES DESIRED 



10 

20 

30 



40 



50 

C— 



NSC=NS*NS 
M H = N S 
M = N2 

CALL CHECK ( E P S , HC , N G, N 0 , IR ET ) 

IF (IRST.E0.1) SSTU3S 
IF (IS IT. EQ- 1) GO TO 20 
CALL READF ( N S . IS A? , E A) 

I? (IDST AE. £q. 0) GO TO 10 
WRITE * (5,1800) 

CALL RDREAL ( A NSR) 

DSTAB= AN SR 
DC 10 1=1, NS 
DESTAB (I) =DST AB 
CONTINUE 
GO TO 30 

CALL SETUP (B A ,G, GAM ,NS , NG, NC) 

CONTIN UE 
WHITE (6,1380) 

DO 40 1=1 , NS 

WRITE (6,1390) (E A (I , J) , J = 1 , NS) 

IF (IDSTAB.EQ . 0) GO TO 50 
WRITE (6,1480) 

WRITE 6,1390) (DEST AB (I ) ,1= 1 , NS) 

CONTINUE 

EIG 2 NS YSTEM OF THE CIEN LOCP DYNAMIC3- 



IF (IOL. EC- 0 . ANC. IQ. EO.O) GO TO 90 
IF 10 L.EQ. 0 . AND. NC. N£. J) GO TC 90 
DO 60 1 = 1 , NS 
DO 60 J = 1 . N S 
60 GN (I, J) =EA (I, J) 

CALL 3 ALA NC { N S , NS , G N .LOW .1 HI GH , D 1 ) 

CALL ORTEES (NS , NS ,LCW, I HIGH , GN, 02 ) 

CALL ORT? AN (NS , NS , LOW, I HIGH , GN. D 2 ,S C) 

CALL BQE 2 (NS r NS , LuW ,IHIGH,GN ,CWR,CVI,SC,IERB) 

IF (IERR.NE.O) CALL EREXIT (NS.GN.I 25 B) 

CALL 3 AL 3 AK ( NS , NS , LOW, I HIGH , E 1 , NS, S C) 

C NORMALIZE AND PRINT OPEN LOOP EIGENSYSTEM 

I W RIT' C ’ = 1 

CALL CNOEH ( C« B , CWI , SC, 51 S , 1 VR ITE , MSQ ,DDD ,D 1 , 02 , WNORB ,W NOBBI , HC , CM, 
1 NO, MS) 

IF (IOL. EC- 2 ) RETURN 

IF (IQ . 2 C.O.OR . (NC. NE.O . OR. IDSTAB.GT . 0 ) ) GO TO 90 
DC 70 1 = 1 , MS 

IF (CWB (I) .LT.O.) GO TO 70 
WRITS ( 5 , 1490 ) 

BETURM 
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70 

80 

90 



100 



110 

120 

130 



140 

150 



160 

170 



180 

190 

200 

210 

220 



230 

240 



250 



260 

270 

280 



290 

300 

310 



CONTI N UE 

IF (IQL • EQ. 3 ) GO TO 130 
DO 80 1=1, NS 
DC 30 J=1,NS 

W11 (1,0) =sc (I.J) 

CALL SlilV (NSQ,W1 1,NS,DDD,D1 ,C2) 
CONTINUE 

IF (IDSTAE.SQ. C) GO TO 130 

?C2fi g * DIAG (DESTAB) * D 

DC 100 J - 1 r N S 
DO 100 1=1. NS 

A A (I, J) =WNORrt ( I , J) * C EST A B (J) 

DO 120 1=1, NS 
DC 120 J = 1 , N S 
DDD=0 . DO 
DC 110 K=1,NS 

DDD=DDD+ A A (I , X) *WiiOKr!I ( K , J) 

DSTORZ ( I , J) = D D fc 
3 A (I , J) = 8 A (I , J) + QCD 
CONTIN UE 

CALL REACH ( N C , NS , IS AH, HO ) 

WHITE (6,1440) 

DC 140 I = 1 , N 0 

WRITE (6,1390) (HC <1 . J) , J = 1 , NS ) 

IF (ia.Ni.1) GC TC 1 So 

CALL i CD E ( W N C H21 , HO , C21, NS , N C , NS , 2 ) 

CONTINUE 

IF (IFDFW.SQ.C) GC TC 170 
CALL 3EACC (NC,NC,D) 

WRITE (6,1470) 

DC 160 1=1. NO 

WRITE (6,1390) (D (I , J) , J= 1 , NC ) 

CONTINUE 
NO 2=0 

IF (NC. EC • GC TC 590 

IF (IOL.FC.Jl GO TC 270 
IF (IR.NE.1.ANE.IE.NE.3) GO TC 210 
IF (ISET.EO. 1 ) GC TC 190 
CALL 3 EACG ( N S , NC , IS AG, G) 

CONTIN UE 

CALL 3FADFB (NC,MS,FEGC) 

WHITE (6,1400) 

DC 190 1=1, NS 

WRITE (6.1390) (G (I , J), J= 1 , NC ) 

IF (ia.Ni.1) GC TC lob 

CALL A CD E ( W NC 581 ,G , BN, NS , NS , NC, 0 ) 

CONTINUE 

GO TO 33 0 

DO 220 1=1, NS 

DO 220 J = 1 , N S 

Rfl (H-MH, J) =0. 0 

CALL H EA CAY (NO, AY) 

DO 240 1= 1, NO 

DO 240 J = 1 , N S 

DDC=0. DO 

DO 230 K= 1 , NO 

CCD=DDD+ AY (I, K) * EC (K,J) 

AA (I,J) = DDD 
WRITE (6,1460) 

DO 250 1= 1, NO 

WRITE (6,1390) (AY (I ,J) , J=1 ,NC) 

DC 260 1=1, NS 
DO 260 J = 1 , N S 
DO 260 K = 1 , NO 

Rfl (I + flH, J)=Rfl (I*fiH, J)+AA (K,I) *BO (K, J ) 
IF (IS ET. EQ. 1) GO TO 280 
CALL 3ZACG ( N S , NC ,1 S AG , G ) 

CONTINUE 

IF (IOL.EC.3) GO TO 290 
CALL REACE (NC,E) 

WRITE (6,1400) 

DO 300 1= 1.NS 

WRITE (6,1390) (G {I . J) , J=1,NC) 

IF (IM.Ni.1) GC TC 3 10 

CALL JCDE (WNORHI ,G ,Ea, NS ,NS , NC,0) 

CONTINUE 



-INV 
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320 
330 
C— - 
340 



350 



360 

C— - 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c— 



IF (IOL. EQ. 3) GO TO 340 
WRITS (6,1410) 

DO 320 Ial.SC 
WHITS (6,13901 (B (I .0) . J= 1 , SC) 
IP (ITF1.SQ.0) GO T 6 3$0 



HO, C M, IP Dr W,D, 38, CC, CP , WR „ W t ,C 
, IT E X ) 



OFEN LO CP TRANSFER FUNCTIONS 
WRITS (6,1500) 

I TF X= 1 

CALL IF (HS,SS,NSC,EA,AA,MC,G,SS,30, 

1 WR,CWI ,SC,JCF ,3ES ,01 .02 , DDD,2E5,ITF1 
IP (IOL. N £. 3) GO TO 260 
IP JNG .EC.O) RETURN 
G C ^0 6 0 0 
CONTINUE 

IF (13. SC. 1. OH. 13. 50. 3) GOTO 500 
CALCULATION OP CONTROL GAI NS; FORMATION OF CONTROL HAMILTONIAN 



T 

I 



-A 



-GM* El *GMT 



-FT 



*** F AND FT ARE THE OPEN LOOP 

DYNAMICS MATRIX AND TRANSPOSE 
***EI IS MCXNC CONTROL WEIGHTING 
MATRIX 

*** A IS THE MSXNS STATE WEIGHTING 
MATRIX 

IS THE NS X HC CONTROL 
DISTRIBUTION MATRIX 



370 

380 
C 



DO 370 1=1, NC 
DO 370 J = 1 , M H 
PRO (1,0) =G (J,I)/B (1,1) 
EC 33 ) 1*1, MH 
DC 330 J = 1 , M H 
HM (I, J *MH) =0. DO 
DO 330 K = 1 , NC 
3 M (I, J ♦ft E) = S' 



c 

390 

C 

400 



410 



420 

430 

440 
C — 



450 
C — 



( I , J ♦ M H ) -G ( I , X) * F?0 (K,J) 

-2NX2N HAMILTONIAN MATRIX 

DIAGONAL ELOCXS Mil AN D M2 2 

DO 390 1-1 , M H 
EO 390 J = 1 , M H 

an ji,j) =sa (i , j) 

B2 (I + MH, J+MH) = -EA (J, I) 

^21 BLOCK 

BM(I+MH,J)=-3M ( I + MH , J) 

M12 BLOCK IS DEFINE! IN LINE 430 AEOVS 

CONTINUE 

I? (IDE3CG.EQ.0) GO TO 410 
WRITE (6,1510) 

CALL 3 A PR NT ( 3 ,M ,2 ,9 .RS, 4 • (9 ( IX. 1PD 13. 6) ) •) 

CALL 3ALANC ( M , M , RM . LOW , I HIGH, El) 

CALL ORTHES (M , 2 , LOW ,IHIGH, RM ,12) 

CALL ORTSAN ( M ,M , LOW ,IH IGH, RM , £2 , X) 

CALL HCB2 ( M , M , L 5 W , I H IG H , RM , W K , W I , X , I ER 3 ) 

IF (IE SR .NE. 0) CALL SREXIT (M , RH, I ERR) 

CALL 3 ALB AK (M ,M ,LOW ,IHIGH, D 1 ,M, X) 

DEBUG DIAGNOSTICS CN EU1EH-LA GRANGE EQU ATIONS 

IF (IDEEUG .EQ. 0) GC TO 430 
WRITE (6,1520) 

DC 420 1=1, M 

WRITE (6,1530) WR (I) , WI ( I ) 

WRITE (6 , 1540) 

CALL R A? ENT ( M , M , M , 9 ,X, 4 , • (9 ( 1 X , 1 PD 1 3 . 6 ) ) 1 ) 

CONTINUE 

IF (IDSTAE .EQ. 1) GC TO 440 
IF (NOe.FC.O) WRITE (6, 1550) 

IF (NOB.NE.O) WRITE (6.1 560) 

IF (NO B. NF. 0) GO TO 750 

CALL 3GAIS ( il , NS , NC , NOB , W R, WI , X , GN , W 1 1 , RM, -21, D1 , CW R ,CWI , S C , M HS , 02 
1 ) 

CHECK EIGENVECTORS 

IF (IDEEUG .EQ. 0) GC TO 450 
WRITE (6.1570) 

CALL RAPSnT (NS, NS, NS, 9, SC, 4, ' (8 ( U , 1PD 1 3 . 6) ) • ) 

CONTIS UE 

RESET FLAG ANE F MATRIX FOR ITERATIVE DESTABILIZATION CASE 

IF (IDSTAE .EC. 3) GC TO 470 
DC 460 1 = 1, NS 
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4 60 SA (1,1) =EA(I, I) -DZSTA3(I) 

13-1 

470 CONTINUE 

C CALCULATION OF FEEDBACK GAIN 

C FEEDBACK GAINS > U = - (BIN VERSE) *G T*GN 

C CALCULATE GT 

DG 490 1=1, NC 
DO 490 J=1,NS 
PRC (I, J) =0. DO 
DG 430 K=1,MH 

480 PHC (I, J) = PRO (I, J) +G { K, I) *GN ( K , J) 

490 FBGC (I , J) =-?30 (I . J) /E (I ,1) 

IF (IDSTAE .EC. i) GC TO 500 

C NORMALIZE AND PRINT OPT. REG. CLOSED LOOP SIG2N3YSTSM 

IWHITE=2 

CALL CNOEM ( CW R , CWI, SC, N S , I SR HE ,NSQ , DD D , D 1 , D2 , WNO R M , W N03MI , F 3GC , 
1AA,MC, NS) 

c THE OPTIMUM FEEDBACK CONTROL GAINS 

500 SHITE (6,1580) 

DC 510 1=1, MC 

510 WRITE (6,1590) (F EGC (I , J ) f J* 1 , NS ) 

C COMPUTE MODAL C MATRIX OPEN LOG? U-INVERS2 SAVED IN WNORMI 

I? (IM . NE. 1 ) GC TO 530 



C IN 
CUE 

C- 



COMPOTING MODAL C RECOMPUTE 0 CPEN LOCP SINCE WNORM USED TO STORE 
U-INV FOR CLOSED LOOP SYSTEMS; SNORM I USED TO SAVE U-INV OPEN LCCP 



DC 520 1=1, NS 
DO 520 J = 1 , N S 

5 20 WNCRM (I, J)=WNCEMI (I, J) 

CALL M I N 7 (NSC / W NCR M , MS . D CD , D 1 , D2 ) 

CALL MODE (W NC EM , £ BGC, A A , NS , N C ,NS , 3) 

53 0 CONTINUE 

C TEE CLOSED LOOP DYNAMICS MATRIX 

DO 550 1=1, NS 
DC 350 J = 1 , N S 
S CM=0 . DO 
DO 540 X=1,NC 

540 SGM=SUM+G (I , K) *F3GC (K,J) 

55 0 ACL (I, J) = EA (I, J) + SUM 

WRITE (6.1600) 

CALL 3 APRNT (MH,MH,MH,5, ACL, 4 , » (5 (1 X , 1 F D 1 3 . 6 ) ) •) 
IF (IR. NE. 1. AND. IR.NE.3) GO TO 590 
DO 5o0 I = 1 , N S 
DO 560 J = 1 , N S 
560 GN (I, J) = ACL (I , J) 

CALL 3ALANC ( N S , NS , GN,LOW ,IHIGH,D1) 

CALL ORTHES ( NS , NS ,LCS , I HIG H , GN , D2) 

CALL ORTBAZ* (NS, NS, LOS, IHIGH,GN,D2,SC) 



CALL H QR 2 (NS , NS , IC V ,IHI G H, G N , CWR ,CW I ,SC, IERR) 
IF (IE RR .NE. 0) CALL EHEXIT (NS , 3N, I3HH) 






CALL 3ALBAK ( N S , NS , LOW , I HIGH , D 1 , NS, S C) 

MALIZE AND PRINT CLOSED LOOP SOBOPT. 



REG. EIGENSYSTEM 



570 

580 
5S0 
6 CO 
610 

620 

630 



— NORMALIZE 
I WHITE = 3 

CALL CNOEM (CS5,CSI , SC, NS ,1 WRITE , NSQ ,DDC,D1, D2,W NORM, W NO RMI, FBGC, 
1 A A , N C ,M S ) 

DO 570 1=1, NS 

IF (CWH]l) .LT.0.0) GOTO 570 
WRITE (5,1610) 

RETURN 
CONTIN UE 

IF (IQ. NE. 1) GO TC 590 
DO 580 1=1, NS 
DO 580 J = 1 , N S 
W11 (I, J) =SC (I , J) 

CALL M IN V (NSC,W11,NS,DDD,D1 ,D^) 

NOB=NO 

IF ( NG . EQ. 0) RETURN 
I? ISET.EQ. 1 ) GO TC 610 
CALL REACG2 ( N S , NG , I G AM , G AM) 

CONTIN OE 

IF (10 L .EQ. 3) GC TO 620 
CALL R EADC (NG ,Q) 

WRITE (6 , 1420) 

DO 630 1=1, NS 

WHITE (6 , 1390) (GAM (I,J) , J= 1 , NG) 
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640 



650 



660 



670 
C — 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c— 



If (IM.Nc.1) GG TC 640 

CALL J1 ODE (WNCRMI,GAM,AA, NS , NS , NG , 1) 

CONTINUE 

IF (10 L .EQ. 3) RETURN 
WRITE (6,1430) 

DO 650 1= 1. MG 

WRITE (6,1390) (Q (I, J) , J= 1,MG) 

IF j( (I Q . EC • 0) .Al^D. ( NG. EQ . 0) ) GC TO 1260 
DO obO I = i , N G 
DO 660 J=1,NS 
FBCjl, J) = 0. DO 
DO 560 K = 1 . N G 

P RC jl , J) = EdO ( I , J) +Q (I,X) *GArt (J,K) 

DO o70 1 = 1 , N S 
DO 670 J = 1 , N S 
CC (I, J) = C.D0 
DO 670 X=1,NG 

CO (I, J) = CQ (I , J) -GAM <I,K) *?RO (K,J) 

IF (IREG .EQ. 1) GO TO 690 

-CALCULATION OF FILTER GAI NS :F 0 R CATION CF ESTIMATION HAMILTONIAN- 



-H0T*3IN*E0 



-Grt*C*GMT 



-FT 



***F AND FT ARE SAME AS FOR 
CCNTFOL HAMILTONIAN 
***Q IS NGXNG STATE DISTURBANCE = 
R COVARIANCE 

***R IS NOXNO MEASUREMENT NOISE = 
aCOVARIANCS 

***BC IS NOXN3 MEASUREMENT MATH I X= 
***GM IS N SX'IG STATE DISTURBANCE = 
DISTRIBUTION MATRIX 



6 SO 
690 
C 



CALL 35AEB (MG ,RC) 
WRITE (6 , 1450 ) 

DO bftO 1=1, NO 



WRITE ( 6 , 1 5 9 0 ) (RC (I,J) , J = 1 ,NC) 
IF (IT F2 .EQ. C) GO TO 700 



-NOISE TRANSFER FUNCTIONS- 

WRIT2 (6,1620) 

ITFX=2 
1 2 E R 0 = 0 

CALL TF (NS, NS, NSQ, ACL, AA, MG, GAM, 3M, NO , HC , CM ,1 ZERO , D , 3B , CC , C? , WR , 

1 WI , CWR ,CiI,SC,JCF.aES,D1 , D2 , DDL, EPS, IT?2,ITFX) 

I? (IREG .aQ. 1) RETURN 
700 CONTINUE 

IF (IREG .EQ. 1) GO TO 930 
IF (IR.LT.2) GC TC 710 
CALL REAEFE (NS,NC,FEG2) 

GO TO 810 
710 CONTINUE 

C THE MEASUREMENT MATRIX (HOT*RI N*HO===>SC 

DO 720 1=1, NO 
DO 720 J=1,MH 

720 PRC (I, J) =60 (I, J) /RC (1,1) 

DO 730 1=1, MH 
DO 730 J = 1 , M H 
R M (1+ 1H , J)=0. DO 
DC 730 £=1,NO 

73 0 HM (I + MH, J) = RM (I+flH, J)-HO (K,I) *ERO (K, J) 

c GM*Q*GMT===>CC 

DO 740 1=1, NS 
DO 740 J = 1 , M S 
RM (I, J) =E A (I , J) 

RM (I + MH, J+MH) = -EA ( J , I) 

74 0 P.rt (I, J + NS)=CQ (I,J) 

GO TO 400 

C-GO 3 ACE TO 45C TO SET UP THE FILTER HAMILTONIAN; CALC. THE FILTER GAINS 
750 CALL RG AIN ( rt , NS , S C , NOB , WR, W I , X , GN ,G M , 3 M , * 1 1 , D1 , CR , Cl , PRO , MHS , D 2) 

c CHECK EIGENVECTORS 

IF (IDS3UG .EQ. 0) GO TO 760 
WRITE (6,1570) 

CALL RAP 5 NT ( NS , NS, NS, 9 , PRO , 4 , • (9 (IX , 1PD1 3.6) ) ») 

760 CONTINUE 

IF (IDSTAE .EC. 1) GC TO 770 

C NORMALIZE AND PRINT OPT. ESTIMATOR El GENS Y ST EM 

I W BIT 2 = 4 

CALL CNORM (C S ,CT , PRC, NS , IW RITE, NSQ, CDD,D1,D2, VNQ Rrt , WNORMI , HO, AA, 



86 



770 

780 



790 



800 

9 10 

320 
C 



83 0 
C“ 



840 

850 



860 

870 



800 

890 

900 



910 



920 
930 
C 



94 0 



950 



1NO,MS) 

DC 780 1=1, MH 
DO 780 0=1, NO 
PEC (I, J) = *60 (J,I) /RC (J,J) 

DO 790 1=1, 3H 
DC 790 J=1,NO 
F3GE (I , J) =O.DO 
DC 79 0 K=1,MH 

FBGE (I , J) =F3GE (1,0) + GN (I , K) * FFC ( K , J) 

IF (IDSTAE .20. 1) GO 70 810 

WRITE (6,1670) 

CALL 3 APRNT ( 3 E , d fi , a H ,5 , G N , 4 , 1 (5 ( 1 X , 1 P D 1 3 . 6) ) * ) 

MBITS (6,1630) 

do 800 1 = 1 , mh 

X (1,1) = DSQRT ( G N (1,1) ) 

WRITE (6,1630) (X |I,I) ,1=1, MH) 

MBITS (6 , 1b30 ) 

DO 820 1=1, MH 

MBITS ( b , 164 0) ( F EG E (I, J) ,J=1,NO) 

COMPUTE MODAL K MATRIX OPEN LOOP U-INV SAVED IN VNORMI 

IF (IN .NF. 1) GC TO 930 

CALL MODE (WNOEMI,FBGE, AA ,KH, SB, NO, 4) 

CONTINUE 

BESET FLAG AND F MATRIX FOR ITERATIVE DESTABILIZATION CASE 

IF (ID ST AB .EQ. 0) GC TO 850 
DC 340 1=1, NS 

DO 840 J = 1 , N S 

BA (I,J) = EA(I,J) -DSTQRE ( I , J) 

15 = 2 

CONTINUE 

DO 670 1=1, NS 

DC 370 J = 1 , N S 

SUM=0. 0 

DC 36) K = 1 , N 0 

SCM=SUM+FSGE (I,K) *8C (K,J) 

FuC(I,J)=£AJI,J) -SUN 
WRITS (0.1639) 

CALL RAPRNT f NS, NS , HS e 5 . PSO, 4 , • (5 (IX , 1 P C 1 3. S) ) •) 

IF (13 . LT. 2) GO TO o9 0 

CALL 3ALANC {NS, NS, PRO, LOW, IHIGH,D1) 

CALL ORTHES ( NS , NS ,LCM, I HIGH , PRO , D2) 

CALL OBT5AN ( NS ,NS , ICS, I HIGH , £SO, D2 , GM) 

CALL HCS2 (No , NS , LOW , Id I Go, PRC,CR,Ct , GM , I ERR) 

IF (IS RP .NS. .G) CALL EBSXIT ( NS , PRO , I S S tt) 

CALL 3 ALEAK (NS, N S , LCW, I 3 IG H , Z 1, NS,G M)' 

WRITE (6.1560) 

NO § KALIS E AND PRINT SUSOPT. ES TIM AT OB SIGENSYSTSK 

I W RITE = 5 

CALL CNOEM (CH,CI ,GM ,NS , I WRIT E , NSQ, D DD, D1 , D2 , W NORM , W NO RMI, HO , A A , 
INC, NS) 

DO 880 1=1. NS 

IF (C3 (I) .LT. 0.0) GO TO 880 
WRITE (5,1660) 

BETURN 
CONTIN OE 
GC TO 900 

IF (IO.2C.0) GC TC 1260 
DO 910 1=1, NO 
DO 910 J = 1 , MH 
PRC (I, J) =0. DO 
DO 910 K= 1 , NO 

PRC (I, J) =1*0 ( I, J) * R C ( I , K ) *P3GE ( J , X) 

DO 920 1=1, MH 
DO 920 J = 1 , M H 
CO (I.J) = G - DO 
DC 9 20 X = 1 , N 0 

CQ (I, J) = CQ (I , J) -FEGS (I, K) ♦PRO (K , J) 

CONTINUE 

THE S MS state AND CONT SOL RESPONSES 

Ifi=IR+ 1 

GO TO (1 CSO, 1 CSC ,S4C ,940) , IR - 

DO 950 1=1, NS 
DO 950 J = 1 , N G 
X (I, J) =0.0 
DO 950 K = 1 , N G 

X (I,J) =X (I,J) +GAM ( I , K) * Q (K,J) 
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960 



970 



980 



990 
1 0 GC 



1010 
1 0 2 C 



1030 

1040 



1 05 C 
1060 



1 07 C 



1080 

1090 

hoc 



1110 

1120 



1130 
1 140 
1150 



1160 



DC 970 1= 1, MS 
DO 970 J=I,NS 
SOfi-O.3 
DO 960 X=1,NG 
S0fl*soa-X (I,K) +GAK (J,X) 
?H0 (I, J) =SU3+CQ (1,0) 

EEC (J, I) = £ EO <I,J) 

CQ (I,J) = SUN 



Q 'JrD *sca 

*21 (I, J) = G3 <I r J) 

¥21 (J,I) =Gtt (J ,1) 

CALL iilNV ( N S C , ¥ 2 1 , N S , D D D , D 1 , C 2) 

CALL SCO 7 (NS ,G2 , ¥21 ,CE ,CI, NS ,Ga,W21 , CR , Cl , PRO , G N ) 
n HITS ( b , 167 0 ) 

CALL RAPRNT ( 2H , 3 H , MH ,5 , G 2! , 4 , 1 (5 ( IX , 1 FD 1 3 • o) ) ■ ) 

¥ RITE (6,1680) 

CC 980 1=1, MH 
X (1,1) =D SORT (G N (I ,1) ) 

WHITE (6 , 1690) (X (I, I) ,1- 1, J1H) 

I? (IQ • EC • 0) GC TC 1260 

DC TOGO 1=1, NC 

DO 1000 J=1,NS 

SUM=0 • 0 

DO 990 K= 1 ,MS 

sua=sua+FEGC ( i, k) *gn (k, j) 

X (I, J) = SUM 
DC 1020 1*1 , MS 
DC 1020 J = 1 , N S 
S UM=0 • 0 

IF (NC- EC-0) GC TC 1020 

do loro k=i,::c 
sua=saa+G (i,x) *x (K, J) 

FRO (I, J) -CQ (1,0) ♦SU£ 

CALL 3 CO V (NS, SC, ¥11 ,CVR , CW I , >1 S , G 23 , W 21 , C B , C I , ? SO , BA) 

I? (NC. EC-0) GC TC 1 C 40 
DC 1030 1=1, NC 
DC 1030 J = 1 , N S 
• 21(1,3) =C. 0 
DC 1030 K= 1 , MS 

W 2 1 (I , J) =¥21 (1,0) ♦ ? 3GC (I , K) *B A ( J , X) 

CC 1060 1=1, NS 
DO 106 0 J = 1 , N S 

sua=o.o 

I? (NC.2C-0) GC TC 1060 
DC 1050 K= 1 , M C 
SUrt=SUa+G (I ,K) *W2 1 (K ,J) 

FRC (I, J) = sus 
DO 1070 1=1, MS 
DO 1070 J=I,MS 

F HC (I , J) = FRO ( I , J) +CQ (I, J) +PRO (J,I) 

PRC (J,I) =FHO (I ,J) 

CALL SCOV (NS, SC, ¥ 1 1,CWR ,CWI,NS,SC,W 11 ,CSR ,CWI ,PRO,CQ) 
DO 1080 1=1, MS 
EO 1080 J = I , MS 
G21 (I,J) =CQ [I , J) 

GM (J,I)=Gft (I, J) 

GO TO 11 CO 
CALL SCOV ( 

IF (NC.EC-t, 

DO 1120 1=1, MS 
DC 1120 J = 1 , M C 
PRC (I, J) = C - D 0 
DO 1110 K= 1 , N S 

PRO (I, J) =PRO (I,J) 4G2 (I,K) *P BGC (J,X) 

CONTINUE 
DO 1140 1=1 , NC 
DC 1140 J=1,MC 
SC (I,J1=C.D0 
EC 1130 K = 1 , MS 

SC (I,J) =SC (I ,J) ♦FEGC (I, K) *PRO (K, J) 

CONTINUE 

IF (IREG -EQ- C) GO TO 1 170 
DC 1160 1=1, NS 
DO 1160 J = 1,NS 
CQ (I,J)=GH(I,J) 

GO TO 11 9C 



*6 A (I,J) - 3A(J,I)+GN(I ,J) 



NS, SC, W 11 ,CWR,CWI,NS,SC,W 11 ,CKR ,CWI ,CQ ,GH) 
* GC TC 1 150 
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1 1 7 C WRITE (6 , 1700) 

CALL 3 A? oST <MK, at, ME, 5, GM, U, ' (5 ( IX, 1?D 13. o) ) 1 ) 

IF (12 .GT.2) GC TC 1 1§0 
CO 1180 1 = 1, 3 a 
DC 1180 J = 1 , 3 H 

1180 CQ (I , J ) *GN (I , J ) +G B (I , J) 

1190 CONTINUE 

WHITE (6,1710) 

CALL 3APoNT ( il E, BE , 3 fi, 5 , CQ, U , ' (5 ( U , 1PD 1 3 . 6) ) • ) 

IF (NC .EC.O) GC IC 1210 
WRITE (6,1720) 

DC 1200 1=1. SC 

1200 WRITS (6,1730) (SC (I ,J) , J=1 , NC) 

1210 DO 1220 1=1, NS 
1 22 C CQ (1,1) =DSQRT (CQ (1,1) ) 

If <&C. EC.O) do TC mo 
DO 1230 1=1, NC 
1230 SC (I, I) = CSQRT (SC(I,I)) 

12UC WRITS (6.17U0) 

DO 1250 1=1, }J5 
IF (I.LS.NC) WRITS (6,1750) 

IF (I.GT.NCj WRITE (6, 1750) 

1250 CONTINUE 
1260 IP (IT ?3 .EQ. C) GO TO 1290 

C FOR3 CCBPENSATOH PEON NEA5 TC INPUT AND COMPUTE T F— • 

DO 1230 1=1, NS 
DO 123 0 J = 1 , N S 
S0M=0. DO 
DO 1270 K=1 , NC 
SUfi=S(J M + FEGc ( I, 5) *HC (K, J) 



c3ii:i)' SC(I ' I) 



1270 

1280 



1290 
C 



1300 

1310 



1 32 C 
1330 



1340 



1350 

1360 



1370 



C 

C570 
13 8 C 
1390 



CQ (I , J ) =ACL (I , J) -SOM 
WHITS (6,1760) 

IT FX=3 
I ZERO= 0 

CALL T F (NS , NS , NSC* CC, A A , NO , FEGE , 3M , NC , FBGC, CM , I ZERO, 3 ,33, CC,CP, 

1 SR , WI, CWR ,CWI ,SC , JCF ,a£S , D1 , C2 , D DD , £ IS , IT F3 , IT ?X ) 

CONTINUE 

COMPETE PSD FUNCTIONS OF THE CO NTBCLLED SYSTEM 

IF (IPS D .EQ. 0) GO TO 13 10 
IF (IYU .IT. 3) GC TC 1300 

CAli PSDCAI (M , NS,R3,X, NC,GW , G 7, FBGC , NO , H Y # H U , HO , ?9 GE, NG , 

1 GAM, ACL ,3A, -5 , WI ,D1 ,02 , JCF , H ES , Q , RC , EE , CC , 1 , IPS D , IN0R3) 

CALL PSD CAL ( 3 , NS ,R 3 , X, N C , GW , GV , F3GC , NO , HY,dY, flO .F3G2, NG , 

1 GAM, ACL f BA, WR , W I ,D 1 ,JCF , RES , Q, RC , S E , CC , 2 , IP SD , IN OR 3) 

GO TO 13 1C 

CALL PSDCAL ( 3 , NS , RM ,X, N C, GW , G V , FBGC , NO , HY , H U , HO , FB GE , NG , 

1 GAM, ACL, BA, WR ,WI ,D 1 ,i)2 , JCF , RES , Q,RC , fi£ , CC, IYU , I PS D , IN0R3) 

IF (ISS . EQ. 0) RETURN 
I? f NC • NE. 0) GO TC 1330 
DO 1320 1 = 1 , NS 
DO 1320 J*1,MS 
ACL (I, J)=BA (I , J) 

CCNTINUl 

CALL MINV (NSQ , ACL ,NS,DDD ,D1 , D2) 

CALL REACW 7 N G , W R ) 

WRITE (6,1770) ( W H (I ) ,1= 1 , NG ) 

WRITE (6 , 17 80 ) 

DO 1340 1=1, NS 

WI (I) =0. 0 

DO 1340 J = 1 , NG 

WI (I) = WI (I) +G AM (I,J) *WR ( J) 

DO 1360 1=1, NS 

C3 (I) =0.0 

DO 1350 J=1 , NS 

CH (I) =CR (iV-ACL (I,J) *WI (J) 

WRITE (6,1390) CR (I) 

DO 1370 1=1, NC 

Cl (I)=0. 0 

DO 1370 J = 1 , N S 

Cl (I) = CI ( I) + FEGC (I, J) *CR (J) 

WRITE (6,1790) (Cl (I) ,1=1 ,NC) 

RETURN 



FORMAT (21, 1P6C14.6,/,2X, 6D 14.6) 
FORMAT (//,5I, 45HCFSN LOOP DYbAMICS 
FORMAT ( 10 (2X , OPE 1 1 . 4)) 



MATRIX ?..,//) 
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1400 


FORMAT 


1410 


FORMAT 


1420 


FORMAT 


1430 


FORMAT 


1440 


FORMAT 


1450 


FORMAT 


1 46 C 


FORMAT 


1470 


FORMAT 


1430 


FORMAT 



//,5X,45HTHE CONTROL DI3TRI3 UTION MATRIX G.. 

//,5X ,45HTHE CONTROL CCST MATRIX ...8.. 

//, 5 X, 45 H PROCESS NOISE DISTRIBUTION MATRIX GAMMA. , 

//,5X , 45HPOWER SPECTRAL DENSITY - PROCESS NOISE 0. . 

// ,51 ,45HMEA3USEMENT SCALING MATRIX H . . 

//, 5X , 45 HPCW ER SPECTRAL DENS ITY -MEASUREMENT NOISE.. H., 

// , 51 , 45 HC IT PUT CCST MATRIX A, 

// ,51 , 4 5 EMEASUP. E M £N T _ F EEDTUR CUG H MATRIX D. 



,//) 

,//) 

,// 

,//) 

,//) 

,//) 

,//) 

f/n 



//,25X.23H DESTABILIZATION CASE.. . . , // , VOX, 39 fiT HE EOLLO 

LUES ' “ ““ 



1490 
1500 
1510 
152 C 

1530 
1 54 C 

155 0 

1 56 C 
1570 
1 58 C 

1 59 G 
1 60 C 
1610 

1620 

1630 

1640 

1650 

1 6 b 0 

1670 
16 6 0 
1 69 G 
1700 
17 10 
1720 
1730 
17 4 C 
1750 
1 76 C 

1770 

1780 
179 C 

1300 



1 WING VALUES WifL EE ADDED DO * N ,/ , 10X , 49 HTHE DIAGONAL OF THE "F" M A 
2TRIX TO DESTABILIZE IT. , / , 1 0 X ,4 1 HOPT IM A L GAINS FOR THE DESTABILIZE 
3D SYSTEM,/, 10X,39EARE THEN USrD AS FIXED SUBOPTIMAL GA INS , / , 1 OX , 28 
4HF0R THE SYSTEM C ALC DLATI CNS . ,//) 

FORMAT (///4 3 H PROGRAM TERMINATING DUE TO UNSTABLE SYSTEM) 

FORMAT <//,2X,3lHCPEN LOO? TRANSFER FUNCTIONS...) 

FORMAT (//,5X,22H EU LER- LAGRANGE SYSTEM MATRIX...,//) 

FORMAT (//, 1X,43HEIGENVALU2S AND EIGENVECTORS OF THE 2N X 2N,/,2X, 

145HEULER-LAGRANGE SYSTEM AFTER 

FORMAT ( IX, 1P2L13.6) 

(IX) 

(//, 2 X , 4 1 B EIG ENS Y S TEM OF OPTIMAL REGULATOR. 



HQR2 . 



,//) 



FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMA] 



( // ", 2 X , 4 1 H El G ENS Y S TE M OF OPTIMAL ESTIMATOR ....IIIII,'//| 
(//,5X,39H EIGENVECTORS FROM RG AIN 



j/'/,1 X',57HTH£ OPTIMAL FEEDBACK GAIN CONTROL 



! ATRIX . 

to unst; 



BLE 



, . F-G*C..,//) 
CLOSED LOOP 



PRIOR TO CNO RM ,//) 

MATRIX. . . C = B INV 

1*GT*S. . . 

FORMAT ( 1 C ( 2 X , 1 P D 11. 4)) 

FORMAT ! // , 2 X f 45 6THE CLOSED LCCP DYNAMICS 
FORMAT ( /// , 6 0 H PROGRAM TERMINATING DUE 

1 SYSTEM) 

FORMAT (// , 21 , 61 H NOISE TRANSFER FUNCIICNS 
loop SYSTEM. . ,//) 

;//,5X, 45HFILTER STEADY STATE 
IX, 2X, 1P6C14.6) 

// r 1 X, 43 HTHE CLOSED LCCP FILTER DYNAMICS MATRIX IS..,//) 

DUE TC UNSTABLE FILTER) 

! E ESTIMATION EERO? P..,//) 

iTIMATION ERROR ,//) 




,3 2H THROUGH THE CLOSED L 



GAINS X . . ,//) 



//, I A, Li.Ui 1U LU'.t' ru>i: 

////,43H PROGRAM TERMINATING D 
//,5X , 45 2 THE COVARIANCE OF T Hi 
// , 5 X , 45 HP MS VALUES OF THE ESI 



<5 (1 X , 1PD13.0) ) ) 

... ... COVA3IANC 



OF THE ESTIMATE XH AT . . , //) 

STATE COVARIANCE MATRIX X=XHAT * P ,//) 

CO NT 50 L CC VARIANCE U=C*XHAT*CT. . . ,//) 



!//,5X , 45 HT H E 
[//, 5X , 45HTHE 
[/ / ( 5 X . 45HTH E 
1P6D14.6) 

f//,2 X , 1 8H STATE RMS RESPONSE, 20X , 20 dCONT ROL RMS RESPONSE,/) 

1 IX, 1PE 15. 7 ,25X, D15. 7) 

r //,5 X ,50'nCCMFENS ATOH TRANSFER FUNCTIONS FROM MEAS. TO INPU 

1T,/,5X,$zH U/Z = -CMSI-F+G*C + K*H)IHV*K. ..,//) 

FORMAT (//, 2X.46HSTEADY D IS TU SEANCE VECTOR fc..,// 

1, 10(1X,1PD12.4/)) 

FORMAT (//,5X, 45HSTEADY STATE VALUES ^OF^ST ATE VAR. ARE ,//) 

, / 



FORMAT (//,5X, 45HSTEADY STATE VALUES OF SI 
FORMAT ( // , 5 X ,47H STEADY STATE CONTROL IS , 
1/ 10 (IX , 1PE12. 4/U 



FORMAT (//,5X , 4^ 
1 , / , 8X , 47 RTO E l 



20 DESTABILIZE 
END 



RENTER THE MAGNITUDE OF THE DESTABILIZATION VECTOR 
ADDED DOWN THE DIAGONAL OF THE "F" -MATRIX, /,3X, 18 HT 
XI- ,//) 
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c== 



10 

20 

30 



SUBROUTINE RAF FNT (N 2AX , 3 ,N , L , A, IDIH , 
REAL** A {NW A X , “ 



DIMENSION 
NU=L 
DC 20 NL=1,N,L 
IF fNU.Gl.N) " 
DC 10 i= i ,;i 
WSITE (6,FMT) 
WHITE (&;30) 
NU=NU+ L 
RETURN 
FORSAT (IX) 

END 



F3T (lllX) 



N 13= li 

(A <I,J) , J = N L, NO) 



FBI) 
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UUUU 



SU3R0UTI NE RG AIN ( M , NS , NC , M C3 , » E , SI , V F , GN , S 11 , TCB , S 21 , LT, C , Cl , CT , M 
1 HS , MT) 

IMPLICIT REAL* 8 <A-H,0-Z) 

PTVP'iqrnV JB 7? / 



DIMENSION 4R (Ml ,SI (flf ,7F (M,3) ,GN (NS, NS) 
DIMENSION Si t (KS ,NS) .TC3 (M, M) ,*21 (NS , NS) , 
DIMENSION C(NS) ,CI (KS) ,CT (NS , KS) 



LT ( NS) ,MT(NS) 



10 



K= 1 
XP = 1 
KN= 1 
N R2E V= 0 
NCPZE7 -0 

I? (K.GT.M) GO TO 210 



CHECK FOR EIGVAL AT CF NEAR J -OMEGA AXIS TO INCLUDE IN E-L EIGSYS 
TORN FIRST ONE POSITIVE AND SECOND ONE NEGATIVE 



20 

30 

40 

50 



60 

70 

C 

80 

90 

100 



C 

1 10 



120 

130 



140 

C 

150 



160 

170 



C 

180 



EIG7R= CAES (SR (K) ) 

IF (EIG7P.GE. 1 .D-10) GO TO 60 

IF (SI (K) ) 4 0, 2C , 40 

NRZEV=N'HZEV+ 1 

IF (N3ZEV.GT. 1) GO TO 30 

SR (K) = EIGVR 

GC TO 80 

SR (K) = - E I G V R 

WRITE (6,290) 

GC TO 1 5 C 

NCPZEV = N CPZ £V + 1 

IF (NCPZEV. GT. 1) GO TO 50 

SR (X) = EIGVR 

SR (K+1) = EIG7R 

GC TO 110 

SR (K) = - E I G 7 R 

SR (X+1) = - E I G V R 



WRITE (6,300) 
GC TO i 8 6 



GC 
IF (SR 
IF (SI 



luc, 70,70 
1 10, 3C, 110 

EIGENVECTOR FOR REAL EIG E NV ALU E , ?OS IT IV E 

IF (NOB. EC. 0) GO TO 100 

DO 90 J=1 ,A 

TCE (J, KP) *7? ( J , K) 

KP=KP* 1 
K— K + 1 
GO TO 10 

EIGENVECTOR FCR COMPLEX EIGENV ALUE, POSITIVE REAL PART- 



IF (NOB. EC. 0) GC TO 130 
DO 120 J=1,M 
FR=VF (J,K) 

FI = -V F (J ,K+1) 

TCB (J, KP) =?3+FI 
TCE (J, X P ♦ 1 ) af S-FI 
KP = K?+ 2 
K=K + 2 
GO TO 10 

IF (SI (K) ) 130 , 150, 180 

EIGENVECTOR FCR REAL EIGENVALUE, NEGATIVE REAL PART 

C (KN) =SR JK) 

Cl (KN) =SI |K) 

I? (NO 3. NE. 0) GO TO 170 

XNS=KN+NS 

DO 160 J=1,M 

TC3 (J , KN S) - VF (Z,K) 

KN = r KN + 1 
K = X+ 1 
GO TO 10 

EIGENVECTOR FCR COMPLEX EIGEN VALUE, NEGATIVE REAL PARI- 
SHES R ( K) 

RI = WI (K) 

C (KN) = HR 
C (KN+1) =RB 
Cl (KN) =BI 
Cl (KN+ 1) =-RI 
IF (NOB.NE.O) 

KNS=KN+NS 
DC 190 J = 1 , M 



GO TO 200 
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I 



T 90 
200 

210 
C 



220 
C 



230 

240 



250 
260 
C 



270 

230 

C — 
2 90 

300 



FE = VF (J,K) 

FI--VF f J # K ♦ 1 ) 

TCE (J,KNS)=FR + FI 
TC3 (J, KNS+1) sFH-?I 
KN=KN*2 
K-K + 2 
GC TO 10 



CONTI N UE 

If (NOB. NE.O) GO TO 



DO 220 1=1, NS 
DO 220 J=1,NS 
W11(I, J) =TC3 ( I , J ♦ N S ) 
CT (I, J) =511 (I , J) 



240 

-FORMATION 



FORM AT 10 N 



OF 511 



Or -2 1 



DC 230 1=1, NS 

DC 230 J = 1 , NS 

5 21 (I, J) =TC3 (I + NS,J+>S) 

IF (NOE. IQ. 3) GO fo 260 

DO 250 1=1, NS 

DO 250 J = 1 , N S 

521 (I, J) =-IC3 (I, J) 

511 (I, J) =TC3 (I+NS ,J) 
CONTIN UE 



INVEST 511 



N SC = N S * N S 

CALL M IN V ( N S C # 5 1 1 , N S , D S TC , L T , MT ) 

CALCULATE THE HG AIN MATRIX 

DC 270 Hal, as 
DO 270 J L = 1 , N S 
G N (IL, JL) =0. DO 
DO 270 KL=1,NS 

GN(IL, JL) =GN(IL,JL) +5 21 ( II, KL ) *5 1 1 (K L , J L) 
IF fNOB.fg.O) SET URN 
DO 280 1=1, NS 
DO 230 J=1 C NS 
CT d,J) =511 (J,I) 

RETURN 



FORMAT (IX, SIR ED IER-UG R ANG E EQU ATI CNS HAVE 
1 1 4 H OR NEAR Zt HO./) 

FORMAT ( 1 X , 4 9 H EUIES-LAG RANGE EQUATIONS HAVE 
1 H EIGEN V AIUES AT 05 NEAR THE J-CMEGA AXIS.) 
END 



A 

A 



REAL EIGENVALUE AT, 
COMPLEX PAIR OF ,40 
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SUBROUTINE MINV ( *SQ , A, N , C, L , E) 
IMPLICIT SEAL * 8 (A-E,0-Z) 

DIMENSION A ( N S C ) , L ( N ) , M ( N ) 

DOUBLE PRECISION A , C , 31 G A , HO L C 
N K = N * N 
C — 1 • 3 DO 
NK=-N 

CO 180 K = 1 , N 
NK=NK+ N 
L { K) =K 

a (K) =:< 

KK=NK+X 
3 IGA =A (KK) 

DC 20 J=K,N 
I2 = N* (J-1) 

DC 20 I=K,N 
IJ=IZ+I 

IT (DAES { EIGA ) -D AES ( A (I J) ) ) 10,20,20 

10 3IGA=A (IJ) 



20 
C — 



30 



70 

C— - 

C 

80 

90 

100 

110 

120 
C — - 



130 

140 

150 
C 



160 
170 
C — - 



L K =1 

a (K i =j 

CONTINUE 



J=I(X) 

I? (J-K) 50 , 5 C ,3 C 

KI=K-N 
DO 40 1=1 ,N 
KI=XI>N 



•I NTSRCKA NGE HCNS- 



40 

r* 


HOLD=- A (KI) 
JI=KI-X+ J 
A (KI) = A (JI) 
A (JI) = HOLE 


50 


I = M(K) 




I? ;i-K) 60 


60 


J P=N* (1-1) 



— INTERCHANGE COLUMNS- 



DO 70 J=1 ,N 
JK= NK + J 
JI=JP+ J 

3CLD=-A(JK) 1 

A (JK) = A ( JI) 

A (JI) = HO L £ 

DIVIDE COLUMN 3Y MINUS PIVOT 

(VALUE OF PIVOT ELEMENT IS CONTAINED IN 3IGA) 

I? (3IGA) 100, 9C, 100 

D=0. 0D0 

RETURN 

DC 120 1=1, N 
I? (T-X) 110, 1 2 C , 1 1 0 

I K=NK* I 



A (IX) =A (IK)/(-3IGA) 
CONTINUE 



DC 150 1=1, N 
IK=NK+ I 
BOLD= A (IK) 

IJ=I-N 

DC 150 J = 1 , N 
IJ=IJ+ N 

IF (I-X) 130, 15C, 130 

IF (J-K) 140, 1 5 C , 140 

KJ=IJ-I+K 

A (IJ) =HCLD* A (KJ) + A (I J) 
CONTINUE 



-REDUCE MATRIX ■ 



KJ=K-N 

DO 170 J = 1 , N 
K J = KJ + N 

I? (J-K) 160, 170,160 

A (KJ) = A (KJ) /BIGA 
CON TIN UE 



-DIVIDE BCV BY PIVOT- 



•PRODUCT OF PIVOTS- 



D=C*3IGA 



C 



A (KK) = (1 • CDO) /EIGA 



•REPLACE PIVOT 3Y RECIPROCAL- 
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CONTINUE 



FINAL HOW ANC COLO JIN INTERCHANGE 



180 
C — 

190 

200 

210 



220 

230 

240 



250 
26 0 



X = N 

K= (K-1) 

I? (K) 260,260,200 

I- L (X) 

IF fl-K) 230,230,210 
JQ=N* (K-1) 

0R = N'*(I-1; 

DO 22d J = 1 ,N 
J K=JO 0 
HOLD* A (JK) 

OI=OR+ J 
A ( 0 X ) = - A (01) 

A (01) = HOLD 

j=« (k) 

IF (O-K) 190 , 190,240 

ki=k-:; 

DO 250 1*1 ,51 
KI-KI>N 
HCLD=A (XI) 

0I=KI- K + 0 
A (XI) =- A (01) 

A (01) = ROLE 
GO TO 190 
K = 9 

RETURN 

SNC 
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10 



20 



30 

40 

50 

60 

70 

30 



90 



100 



110 

120 

130 



140 

150 

160 

170 



180 

190 




SUEHOOTINE SCCV (NL, «L,NLI, VL 1 ,VL2, N R , W B , W RI , V R 1 , VR2, 

REAL*** 7L 1 (NL) ,712 (ML) , WL (NL, NL) , '^LI (NL,NL) , X ( N L , N R ) , * nn/ 

1 VR 1 (NR) ,VR2 (MS) ,'ifi (NR ,NR) ,M FI (NR, NR) 

HEALTH A,e,C, D,X1,K2,K3,K4 
DO 20 1=1 , NL 
DO 20 J = 1 , N R 
X (I,J) =0. 

DO 20 11*1, NL 

X (I,J) =XJI,J) + SLI(I,II)*Q (II , J) 

DC 40 1= 1, NL 
DC 40 J= 1 ,NR 
Q (1,0) =0. 

DO 30 JJs1,M3 

G (I,J) =Q (I,J) + X <I,JJ) *SRI (J, JJ) 

CO NT IN Oc 
1=1 

IF (VL2 (I) ) 60,110,60 
J= 1 

IF (VR2(J)) 3C ,90,80 
A= VL1 (I) +VR 1 (J) 

B=-2.* VL2 (I) * Vfi2 (u) 

C=A**2 +VL2 (I) * * * + 7 H 2 (J) **2 
D=C**2-3**2 
X 1 = A*C/D 

K2=- ( V 32 (0) *C + VL2 (I) *B) /D 
K3 = - (7 32 (J) * E 4 VL2 (I) *C) /D 
K4*-A*B/D 
11=1+1 
J1=J+1 
X (I, 

X ■ 

X i 

X ( J. I ,U 

J =0 + 1 
GO TO IOC 
A=VR1 ( J) +VL1 (I) 

B=A**2 +VI2 (I) * *2 
K 1= A/B 
K 2 = VL2 (I) /3 

X (I,J) = K1*Q(I,J)-X2*Q(I+1,J) 

X (1+ 1 , J) = K2*Q (I, J) + K1*Q (1 + 1 , J) 

J=J+1 

IF (J. LE • NR) GC TC 70 
1 = 1 + 2 
GO TO 160 
J = 1 

IF (VR 2 ( J) ) 1 3 C , 140, 130 

A = VP 1 fj) + VL1 (I) 

E=a**2+VR2(J) **2 
K1=A/9 
K2=VR2 (J) /B 

X (I,J) =K 1*0 (I , J) -*t<2*C (I , J + 1 ) 

X (I, J+ 1) =K2*Q ( I, J) +K 1*Q (1,0 + 1) 

J=J>2 
GC TO 150 

X (I , J) =Q (I,J) / (VR 1 (J)+VL1 (I) ) 

J=J+1 

IF (J.L2.NP) GO TC 120 
1 = 1+1 

IF (I. LE.NL) GC TC 50 
DC 170 1=1, NL 
DO 170 J=1, NR 
Q (1,0) -0 . 

DO 170 11=1 , N I 

Q (I,J) =Q <I,J) +«L(I,II)*X (II,J) 

DO 190 1=1, NL 
DC 190 J*1,NR 
X <I,J) =0- 
DC 180 J J= 1 , N 3 

X (I,J) =X (I,J) +c (I , J J) *W R ( J, J J ) 

CONTINUE 
RETURN 
END 
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ii me 



nnonnnnnnnnnn 



C== 



SUBROUTINE MODE ( WNC 5H, G , GNOR M, NS , Ml ,N2,ICCM) 



WNORM TR ANSFORMATICN MATRIX U C5 0-INV 

1.2 NO . OB STATE 

NC NO. OF INPUTS OP OUTPUTS 

ICON C0NT3CI FLAG TC INDICATE WHICH TR ANS FOB 21 AT ION 

0 = MODAL G 

1 = MODAL GAMMA 

2 = MODAL H 

3 = MODAL C 

4 = MODAL K 

5 = CONTROL EIGENVECTOR MATRIX 

6 = MEASUREMENT EIGENVECTOR MATRIX 



10 

20 

30 

40 

50 

60 

7 3 

30 

90 



100 

110 

120 

130 

140 

150 

160 



C 

170 

180 

190 

200 

210 

220 

230 

240 



IMPLICIT HEAL *8 (A-B.C-Z) 

DIMENSION WNORH (SS ,NS) , G < SI , N2) , GNOR K ( N 1 , N2) 
CO 10 1*1, SI 
DO 10 J= 1 ,N2 
GNORM (I, J)=0. 

I POINT = ICCN * 1 

GO TO (2G,20, 3G, 90,20,90 ,90) , IPCINT 
DC 30 J=1,N2 
DO 30 1 = 1 ,MS 
DC 30 K= 1 *js 

Gsoa; j . (i , j) =g:icbh ji, j) -t-aHOEa (i , ?) *G(K , j) 

GC lO (40,70, SC, 3C, 80) , IEOIHT 
WRITS (6 , 170) 



« r. x x _j \ ^ \ 

CO 60 1=1 .NS 
WHITE {6 , 230) 

RETURN 

WRITS (6,180) 

GO TO 50 
WHITE (6,240) 

GO TO 50 
DC 100 J = 1 , N S 
DO 100 1=1, Ml 
DC ICO 5-1, NS 
GNORM (I. J)=GNCRM (I, J 
GC TO ' (‘ ‘ ‘ K ^ 

WRITE ( 

TO T 



(GNCRH (I, J) ,0=1 , M2) 



GC TO T 5 d 
WHITE (6 .200) 
GO TO 1 5 C 
WHITE (6,210) 
GO TO 15 C 
WHITE (6,220) 
DO 160 1=1. Ml 
WHITE (6,230) 
RETURN 



. II, J)+G (I .X) *WNCRM (K , J) 

1 10,1 10, 11 C, 1*0, T 10, 13 C, 14 0) , IPOINT 
6 , 190) 



(G NORM (I, J) ,J= 1 , NS) 



FORMAT 

FCBMAC 

1 ?6hmat 

FORMAT 

FORMAT 

FORMAT 

FOHMAT 

FORMAT 

ENC 



(//, 5X,45H MODAL CONTROL DISTRIBUTION MATRIX 

(// , 5 X , 5CHMOCAL PROCESS NOISE DISTRIBUTION MATRIX 



TI*G- . ,//) 
. . .TI*GAC. 



(//, 5X , 45HMCDAL MEASUREMENT SCALING MATRIX 

//,5X,45HTHE MODAL CONTROL GAINS 

(//, 5 X , 45 HOC NTRO L EIGENVECTOR MATRIX 

(//,5X,45HMEA2UR£MENT EIGENVECTOR MATRIX. • 
{ IX, j2X, 1P6D14.6) ) 

//,DX,45b'MCD AL FILTER STEADY STATE GAINS. 



H { 3AR ) *T. 

.. ...C*T. 



C*M. 

H (3AR) *H. 



TI*K. 



,// 

,/A 

,//) 

,//) 

,//) 
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C== 



SUEROUTINE CMC EM 
1C,CM,N 1, N 2) 



(WZ ,WY , VEC,NS,IWRITE,NSQ,DDD, D1 ,D2,WN0RM 



wnorki,h 



*Z (I) 
WT (I) 
VEC 
NS 



SEAL PART OF I-Tfi EIGENVALUE 

CCS FLEX PART OF I-TH EIGENVALUE = 

MATRIX 0? RIGHT EIGENVECTORS STORED IN REAL FORE 
FF.CM HQR2 = 

NO. OF STATES 



IW SITE FLAG TC CONTROL FORMATS FOR DIFFERENT EIGHENS YSTEMS= 



WNORM NORMALIZED MATRIX U CF RIGHT EIGENVECTORS STORED 

BY COIUMNS IN FEAL FORM 

WNORMI O-INVERSE 2*C0NGUGATE Or LEFT EIGENVECTORS 
STORED BY ROW IN REAL FORM 
NSQ, CCD, D1 ,C2 - ARGUMENTS PASSED TO MINV 



10 

20 



30 



IMPLICIT REAL * 8 (A-H,Q-Z) 

SEAL* 3 FIELD , CCMM A , S FMCOL ,RIGHI,FMT 
DIMENSION WZ (NS) , 5Y (NS) , V EC ( NS, NS) .1 NORM (NS, NS) 
1E(6),D1(NS),D2jNS),FMT(14) c HO (N 1 , N2) , CM (N 1 4 N2) 
DATA FIELC/5HE12.5/,C0MMA/^H, 1 , 1 , /, S EMC0L/5H , * : 
16H(1X,1P,13*1H / , S E M E ND / 4 H , 1 : »/ 

NORMALIZE COMPLEX EIGENVECTORS BY LARGEST 

KK = 0 
L H=0 
LC = 0 

DC 50 K= 1 , NS 

IF (XX.2Q.1) GC TC 40 

IF (DABS ( WT ( K) ) . LI. 1.D- 10) GO TO 50 

LC=LC* 1 

EM AX=0 . DO 

DO 20 1=1, as 

CMCD=7 EC (1,5) * * 2 ♦ VEC (I, K+ 1) ** 2 
IF (CMOD-EMAX) 20,10,10 
EM AX=CMOC 

a=i 

CONTINUE 

V tt R= V S C (M,X) 

7MI= VEC (M ,X+ 1 ) 

DO 30 1=1, NS 

V R = VEC (I ,R) 

V 1= V EC (I,X+1) 

VECRN= (VR*VMR+VI*VMI) /E MAX 

V ECIN= (-VR*7MI+VI*VMR)/EM AX 
WNORM (I , X) =VEC5N 

WNORM ( I , K ♦ 1 ) = VECIN 
CONTINUE 



, W NO SMI (NS, NS) , STOP. 
• , / , SIS HT/1 H) / , FMT/ 
ELEMENT 



K K = 1 

GO TO 50 
40 K X =0 

50 CONTINUE 

c NORMALIZE REAL EIGENVECTORS BY THE TOTAL LENGTH 

DC 30 K = 1 , NS 

IF (DABS (WY(K) ) .GE. 1 .D-10) GO TO 30 
L R=LR> 1 



60 



90 

100 



1 10 



120 



R EMOD = 0 .DC 

DC 60 1=1, NS 

REMOD= VEC (I , X ) **2 ♦ RE MOD 

BMOD= DSQ5T (REMOC) 

DO 70 1=1, NS 
R V EC= V 2C (1,X) /EKCD 
VNCRM (I, K) =RVEC 
CONTI NUE 

CONTINUE — 

GO TO (9C f 10G,110,12C,130), IWFITE 
WRITE (6.320) 

GO TO 140 
WHITE (6,330) 

GO TO 1 4 C 
WRITE (6.340) 

GO TO 14 C 
WRITE (6,350) 
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130 

140 



C 



150 



160 



170 

190 



190 

200 

210 

220 
23 G 
240 

250 

260 

270 

280 

C 

290 

300 



310 



GC TO 140 
WRITE (6 ,360) 

KK=0 
NPSTW^O 
NFS TV- 1 
DO 180 1=1. NS 
IF (KK.EC.l) GC TC 170 
IF (DABS (WY (I)) .GT. 1 .C- 1 0) KK- 1 
PRINT OUT NO MOnE TRAN 6 WORDS, NOT SEPARATING COMPLEX 
IF (N? RTW.LT. 5 .OR. (NFRTW. £Q. 5 . AND. KK . EC . 0 ) ) GO TO 150 
FMT (NFMTW+1) = R IGHT 

WRITE j(6,FMT) (STCRE (J) ,J=1,NFFTW) 

NPRTW-0 

nfmtw= 1 

NPRTV=NPRTW+ 1 
Nr MTV = NF STW ♦ 1 
IF (KK.EC.l) GC TC 160 
STCRE (NPFTV) = WZ (I) 

FMT (MFMTW) =FIZID 
NFMTW= NFMTW+ 1 
FMT (NF MT 7) =3E MCCL 
GC TO 180 
STORE (NP HTV) = W2 (I) 

FMT (Nr MT V) =FIE1C 
FMT (NF ST W+1) =CCMMA 
STORE ( NP RTV ♦ 1 ) = WY (I) 

FMT (NF MT V + 2) =FIELD 
FMT (NFMTW+3) =SEMCCL 
N F MT W= NP MTW+ 3 
N P3T V = NPBTW+ 1 
GO TO 1 8 C 
K K - 0 



IX, 1 P D 1 3 . 6 ) ) ») 



CONTIN US 

FMT (NFMTW) =SZMENC 
FMT fNFMTW + 1) =RIGHT 

WHITE (6,fMT) fSTCRE (J) ,J=1,NF?TV) 

IF (I VP.I TE . N E . • ) EG TO 1 c: > 0 
WRITS (6,370) 

GC TO 2 0 C 
WHITE (6 ,330) 

CALL RAPBNT ( N S , NS , N S .6 , S NO R M , 4 , 1 (6 ( 

GC TO (230, 21 C , ‘ 10, 2 20, 220) , IVRITS 
CALL MODE ( W NORM ,HO , CM, NS , N 1 , N2 , 5) 

GO TO 23 C 

CALL M CD E ( W N C SM , HO , CM , N S , N 1 , S 2 , 6) 

GO TO (24C, 250, 26C ,270,280) , IWRIT2 
WRITS (6,3<50) 

GO TO 290 
WHITS (6,400) 

GO TO 29 C 
WHITS (6,410) 

GC TO 290 
WHITE (6 ,420) 

GO TO 290 
WHITE (6,430) 

SAVE U-INVERSE CPEN LOOP IN WNORMI 

IF (IV3ITF .GT. 1) GC TO 310 
DO 300 1=1, NS 
DC 300 J = 1 , N S 
WNORMI (1,3) = W N C F M <1 ,3) 

CALL M I N V (NSQ, WNORMI, NS ,DDD , D1,D2) 

CALL R A ? 3NT ( N S , N S , NS ,6 , W NO RHI ,4 , • (6 ( 1 X , 1 P D1 3. 6) ) ' ) 
BETU2N 

CALL MINV (NSC .SNCEfl ,NS , DDD, D1.D2) 

CALL R APB NT ( N S , NS , N S ,6 , WNO RM , 4 , • ( 6 ( 1 X , 1 PD 1 3 . 6 ) ) 1 ) 
RETURN 



EIG7 AL 



C 

320 

330 

340 

350 

360 

370 

380 

390 

400 



FORMAT 

FC3MAT 

FORMAT 

FOBMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 



( / / 1 X , 4 2RO FEN LOO? EIGENVALUES D ET ( S I- F) . . , //) 

(//IX, 46HC-LCOP OPTIMAL REG. E-VALUES 3 3T ( S I-F+G* C) ..,//) 

(//IX, 46HC-LOCP SU30PT. REG. E-VALUES DET (31- F+G*C) . . ,//) 

1 // 1 X , 4 6 HC-LCCP OPTIMAL EST. 2-V ALU ES . . . D ET ( SI- F + K * H) ..,//) 
//IX, 46HC-LOCP SUBOPT. EST. E- VALUES. .. 3ST (SI- ?♦ K*H) .. ,//) 

//I X , 4 6HC PE N LOOP RIGHT EIGENVECTOR MATRIX T ,//) 

//IX, 46BC-LOO? RIGHT EIGENVECTOR MATRIX M ,//) 

// 1 X , 4 6HCES N LOO? LEFT EIGENVECTOR MATRIX T-INV..,//) 

//IX , 4 6 HC-LCCP OPT. REG. LEFT E-VECTOR MATRIX. . M - 1 N V . . ,//) 
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410 

420 

430 



FORMAT 

FORMAT 

FORMAT 

ENC 



(//IX, 46HC-LOO? SUBOPT-REG. LEFT E-VECTOR MATRIX. . « 
i//\ X, 46HC-LCCP OPT. FILTER L EFT E-VECTOR MATRIX. . M 
(//IX , 5 1HC-LGO? SU50PT. FILTER LEFT E-VECTOR MATRIX 



INV,//) 
INY,//) 
. M-IRV. 
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c== 



SUE30UTI HE TP ( S , S f! , NSQ . A , A A , M ,3 , B« , l , C , CM , I FD PM , D , BB 
1 S va ,EVI,?H , ‘l,3C,JCF ,3ES ,D1 ,32 , EDO, EPS, XT F, ITFX) 

IMPLICIT 9EaL*S (A-a,C-?,) 

DIMENSION A (S , S) , A A (N-N) , B ( N , E) ,3(1(3 ,3) ,C (L, H) ,CM (L, »| ,D (L ,3} ,3B <N 
,CC(M) ,CE (N) , EV3 |H) , EV I (N) , PS |N) ,?I (}!) ,SC (N ,N) , JCF ( N) ,SES (!l) ,D1 (H 



F 33 r CC , C ? 0 



1) ,CC(N) 

save' comp It 
if <itf 



c— 

10 

20 



30 



40 

50 



60 



TAIION 
(ITFX .2Q. 



CN CL AND CL SYS WITH MODAL WORK DONS IN OPTSTS- 
GO TO 50 



I? |IT?X .2Q • 2\ GO TO 10 

CALL ?CL£S ( N , N« , A, AA,M , 3,L, C f ?R,?I, D1 , D2 r JCF,5C) 

COMPUTE ttCCAL i! Aral CIS FOR RESIDUES 

DO 20 1=1,11 
DO 20 J = 1 , tf 
AA (I,J)=SC(I, J) 

DO 30 1=1 ,L 
DC 30 J=1,:J 
CM (I,J) = C.D0 
DC 30 K=1,N 

CM (I, J) =C3 (I, 0) *C |I , K) *AA IK, J) 

CALL 3INV (NSQ,AA,N,CDD,D1,D2) 

DO 40 1=1, N 
DO 40 J=1 ,3 
E3 (I, J) =O.DO 
DO 40 K= 1 , N 

3H (I,J) =Efl (I, J) >AA(I,K) *B (K, J) 

CONTINUE 
DO 60 1=1,3 
DO 60 J = 1 ,L 

IF (IT F • NS. 3) CALL ZEROS ( I , J , IFDF N , SI , MM , A , A A, £1 , B , L, C, D , BH , CC , CP 
1 ,2Va,2YI ,D1 , D 2 ,SPS) 

IF (IT? . NS. 2) CALL RES ID ( I , J , N , JC F , 3 , 3 3 , L ,C3 , P R , PI , RES , 3B , CC , 1) 

CCNTIN UE 

RETURN 

END 
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c 




DO 10 J=1,N 

10 AA (I, J ) = A (I, J) 

CALL BALANC ( N 3 , N , A A ,LO W - III IG F , D 1 ) 

CALL OBTHES (NM,N, LOW, THIGH, A A , D 2 ) 

CALL OHTRAN ( N Z , S ,LO '* , I HI 3H , A A ,D2 ,SC ) 

CALL HCH2 (N:i,N,LOW, 1HIGH,AA,EVR,2VI,SC,ISSR) 

IF (I2P.R -NS . C) GC TO 30 

CALL 3 ALB AK ( N fl , N r LO * , I H IGH , D 1 , N , SC) 

WHITS ( b .40) 

CC 20 1=1, N 

20 WRITS (6,*Q) E7R (I) , EVI (I) 

HEiaP.N 

30 WHITS (5,60) 

5ZTUBN 



C 



40 FCB3AT 
50 FORMAT 
60 FCHHAT 
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S CEEOUTINZ ZEROS (K1 , K2 , IFD FW , N , N N, A , A A , N, 3 , L, C , D , B d, CC, C? , E V3 , 271 
1, D 1 , D2, ZP S ) 

I HPLICI T REA L * 8 ( A -H , C-Z) 

DIHENS ION A <N , N) , A A (N,N) , B ( N , N) ,C(L, N) , D ( L r Cl ) , B3 (N) , CC (N) , CP ( N ) , E7 
1R (N) ,2 71 ( N) ,01 (M ,02 (Ml 
DOUBLE PRECISION SCL ,DA3S 
DC 10 1=1, N 
32 (I) =3 (I x X1) 

CC (I) =C (Xl f I) 

DC 10 J=1,J 
AA (I, J) = A (I, J) 



10 



20 



30 



40 

50 



60 

70 

80 

C 

90 

100 

110 

120 

130 

140 



WRITE (6, 90) XI, K2 
IF (IFDFW .Eg. 0) GO TO 20 
H = D (X2 ,K 1) 

IE (DABS (n) . Lc. EPS) GO TO 20 
J J = N 

GC TO 50 
NN=N- 1 

DC 30 1=1, NN 
H = SCL ( N , 22, CC) 

CALL CCOHP ( N , N H , A A , CC, CP) 

IF (DABS (H) .GT .EPS) GO TO 40 

CONTINUE 

H = SCL (N ,£E,CC) 

WHITE (o , 100) H 
GO TO 70 
JJ = S T -X 

WRITE (6,110) u I , H 

CALL A CO HP (N,NK,AA,E3,CC,H) 

CALL 3ALANC { N H , N , A A , LOW IH IG F , D 1 ) 

CALL OBT blS ( N H , N ,LO W , IHIGH , A A , C2) 

CALL 30 R (NN , N , LC 5 .1 HIGH , AA , E V 5 , EVI , I ES F) 
IF ( IE 22 . N E . 5) GO TO 80 
W R I " 



TE (6.120) 
50 1*1,3 



DC 

WRITE (6,130) £7R |I) ,EVI (I) 

RETURN 

WRITE (5,140) 

RETURN 



IH:) 



FOSNAT {///, 1 7 H TF FOR INPUT NC.,I3„15H AND OUTPUT NO. ,13, 

FORHA7 // c oX,27HNG FINITE ZEFCS. TF GAIN =,E12.4) 

F03HAT (/, jX, fOHCSDER OF NUMERATOR =,I3,0X,99T? GAIN =,S12.4) 
FORMAT (/ ,3X, 57HNUH2SATOR EIGENVALUES (INCLUDING EXTRANEOUS ZERO 7 
1 ALU SS) : ) 

FORMAT ( / ,4 X , 1H(,F13.6,4H)+J(,F13.6, IH) ) 

FOR NAT ( 52H FA ILU32 IN HQR CALCULATING TRANSFER FUNCTION ZEROES) 
END 
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owe/* u 

^hUOOHNG II 
CJI3M UU-rti II 



G3TI32 AC02!? (S,SS,A f a,C f t) 
A . 2 . C . r. 

NS ION A (5 2,3) f 5 { !J) f C ( N 1 ) 

0 1=1 ,N 

o j= i r :x 

J) =A (I,J) -2 (I) *C | J) /H 
2N 



104 



SUBROUTINE CCCMP (N,N»,A, 

R S AL *8 A.C.CC 

disess ion a (::a ,N) ,c (N) , cc (N) 

DO 10 1=1, N 
CC (I) =0. 

DO 10 J= 1 ,N 

10 CC (I) =CC (I) *C (J) *5 (J,I) 

DO 20 1=1, N 
20 C(I)=CC(I) 

RETURN 

END 
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nil 



awnonHwc n 



NOTION SCL ( N , 2 ,C) 
AL*$ 3 , C r SCL 
jE NSIO N 2 (N ) ,C (N) 

10* 1=1 ,N 
L=SCL + C (I) * S (I) 
TU3N 
E 
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SUBROUTINE RE SID (K 1 ,X2 , N , JC? , fi, 3M, L , CM , ?S, ?I, a BS , 3B,CC, IPT) 
IMPLICIT BEALES (A-H,C-Z) 

i DIMENSION JC? (N) , EM ( N, M) , CM ( L , N) , ?R ( N) , PI ( N) ,333 (N) ,33(N) ,CC(N) , P3 

cIta s n/8e*sin (b*t/, a i/bh */, a 2/sazzp (A*T) / f sd/i h) / 



10 

c— 



20 
c — 



30 



40 

50 

c— 



c— 

60 



70 

80 

C 



90 



DATA Z£RC/0.D0/,T1/4H*T**/, B L A NX/8H /,C5/3H*CGS (B*T/ 

TEMPORARY mod till jcf is calculated 

DO 10 1*1, N 
JCF (I) =0 

TEMPORARY M CD 

IF (IPT .10. 1) W FIT 2 (6,170) 

DO 20 1=1, N 
BE(I) = EM ( I , K 1 ) 

CC (I) = CM (K2,I) 

loop THROUGH ihe poles 

1=0 

1=1*1 

IF (I .31. N) GO TO 160 

IF (JCF(I) . EQ. 1) GO TO 60 

IF (DA BS ( PI ( I ) ) .LT. 1.0-10) GC TO 50 

COMPUTE SIMPLE COMPLEX POLE RESIDUES AMD PRINT BOTH 

RES (I) =CC(I) *BE (I ) + CC (I * 1 ) * B B (1+1) 

3ES (1+1) =CC(I) *E£ (1+ 1)-CC (I + 1 ) *3B (I) 

IF (IPT .30. 0) GC TC 40 

PRT (1) =BLANK 
PBT (2) = R 2 

IF (PI (I) .EQ. 0 . DO ) ?2T (2) =3IANK 
FRT (3) =CS 
PRT (4) =ED 

WRITE (6 , 180) PR (I) , PI (I) ,RSS (I) , (PRT (J) , J = 1 ,4) 

1 = 1+1 
PBT (3) =SN 

WRITE (6,130) PR (I) , PI (I) ,RSS (I) , (PRT (J) ,J=1 ,4) 

GC TO 30 
1 = 1+1 
GC TO 30 
CONTIN UE 

COMPUTE SIMPLE REAL POLE RESIDUE 



PR (I) ,?I(I) ,RES (I) , (PHT(J) ,3=1,4) 



RES (I) =CC (I) *3E (I) 

IF (IPT .EQ. 0) GC TO 30 
PRT (1) =R 1 
PRT 2 =R 2 
PBT (3) = BLANK 
PRT (4) =3 LANK 
WHITE (6,130) 

GC TO 30 

LOOK AHEAD TC DETERMINE SIZE OF THE JORDAN BLOCK 

K= 1 

KT= N- 1 

DO 70 J=I,KT 

IF (JCr(J) .EQ. 0) GOTO 80 

K = K + 1 

CONTINUE 

IF (DABS (PI (I)) .IT. 1.D-10) GC TO 110 

COMPUTE REPEATED COMPLEX PCLE AND PRINT OUT ALL FOUR- 



K = 1 

RES (I) =CC (I) *EE (I) +CC(I+ 1) *B3 (1+ 1) +C C (1 + 2) *3 3 (1+2) +CC (1+ 3) *3 8 (1+ 1) 
RES (1+ 1) =CC (I) *BB (1+ 1) -CC (1+ 1 ) *33 (I) +CC (1 + 2) *B8(I+3)-CC(I+3)*6B(I* 
12 ) 

RES (1+2) =CC 
RES (1+3) =CC 
IF (IPT .EQ 
PRT 
PBT 



(I) *33 (I+3)+CC (1+1) *BB (I* 2) 
(I^*B3jI + 3)-CC (I + 1)*B3(I+2) 



)) GC TO 100 



(IPT . 

1) =R 1 
(2) =H2 

IF (DABS (PR (I) ) .GT. 1.D-10) GC TO 9 0 
PBT (1) = BLANK 



PBT (2 =8 LANK 
PBT (3) =CS 
PRT (4) *ED 
W HIT E (6,180) 
PRT (3) = S N 
1 = 1+1 

WRITE (6.180) 
PBT (1) =T 1 
PRT (2) =R 2 



PR (I) ,PI(I) ,RSS (I) , (PRT(J) , J=1 ,4) 
PH (I) , PI (I) , RES (I) , (PRT ( J) , J=1 ,4) 



107 



100 

110 



120 

130 



IF (DABS (£3(1) ) .LT. 1. D- 10) PRT ( 2) = BLA NK 
PET (3) =CS 
1 = 1+1 

WP.ITS (6,190) PH (I) ,PI (I) ,RSS (I) ,PBT <1) ,K, (PHT (J) ,J = 2,4) 
PET (3) =S N 
1 = 1+1 

WHITE (6,190) PS (I) , PI (I) ,RES (I) , PRT (1) ,K, (PRT (J) , J=2, 4) 
GC TO 30 
1=1+3 
GC TO 30 

-COMPUTE REPEATED HEAL POLE RESIDUE AMD PRINT OUT ALL K OF 
CONTINUE 
KT=I+K-1 

;in=o 

DO 130 J=I,KT 
N N=NN + 1 
RES (J) =ZZSO 
DG 120 J J = J , K T 

RES (J) =RES (J) + EB ( 0 J ) *CC (JJ-NN + 1) 

CONTINUE 

IE (I? T .EQ. 0) GC TC 150 
N N=0 

PRT (1) =T 1 
PHT (2) =R 2 



* a EM 



PRT (3) =B LANS 
FST (4) = EL ANK 
DO 140 J=I.ST 



140 

150 

160 
f* — 


WRITE (6 
N N = NN+ T 
GC TO 30 
I = KT 

GO TO 30 
COHTIM UE 
RETURN 


170 


F CBM AT' (, 


180 


IE 3 I D 
FORMAT (, 


190 


FORMAT (, 




1 2 AS, A 1 ) 
ENC 



1 = 2 , 4 ) 



' ,19 .7HP.E AL (A) , 126 -7 HI HAG (5) ) 



1 5 H R 

A IV 

^ t 2 X t 
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C — 3 



c — 
10 



20 



30 
40 
C — 
50 



60 



70 



90 



C 

9 0 
100 



110 



120 
C 

130 
C — - 
140 



S U EHO'JTI NE 3 ALANC (N H, M , A r LOW , IGH ,SC ALE) 
I N TEG 2 3 I , J , K , L , 3 . N , 0 J , N a ,IG H , LC W , I E XC 
REALM A(NH,N) ,SCALMN) 

REALM C,f,G,R,S,E 2 , RADI X 
REALM 0i£3 
r o/? tp at wrrn v 

DATA aADIX /242 1 C 0000 C 0000000 / 



3 2 = RAD IX * RADI X 

K= 1 

L=N 

GO TO 60 

IN-LINE PROCEDURE FOR ROW AN C COLUMN EXCHANGE 

SCALE (H) = J 

IF (J .EC. *) GO TO 40 
DO 20 1 = 1 ,L 
F= A (I r J) 

A (I,J) =A ( 1 , 3 ) 

A (I, a) = f 
CONTINUE 
DC 30 I=K,N 
F = A (J I) 

A (J r lf = A (3,1) 

CCNfllfUE 

GO TO ( 5 C. 90 ) , IFXC 

’ SEARCH ?C 5 30 * S ISOLATING AH EIGENVALUE AND PUSH THEM DOWN* 
GO TO 230 



.EC. 1) 



GC TO 70 
. NE . 0.0D0) GO TO 30 



IF (L 
L= L- 1 

DO 60 J J = 1 , L 
J=L+ 1 -JJ 
DO 70 I = 1 ,L 
IF (I .EC. J) 

I? (A(J t I) 

CC NT IN U i 
fl = L 

I EXC= 1 
GO TO 10 
CONTINUE 
GO TO 100 

SEARCH FOR CCLU 3 NS ISOLATING AN EIGENVALUE AND PUSH THE a LEFT- 

K = K+ 1 

DC 120 J=K,L 
DO 110 I =K , L 
IF (I .EC* U) GC TO 110 

■ " , N E . O.ODC) GO TO 120 



ufe J) 



IF i 1 

CONTINU 
a=K 

I EXC= 2 
GO TO 10 
CONTINUE 

NOW BALANCE THE SUE 3 ATEIX IN ROWS K TO L- 

DC 130 I = K , L 
SCALE (I) = 1.000 

ITERATIVE LOOP FOR NORil REDUCTION 

NOCONV=. FALSE. 

DO 220 I = K,L 
C = 0 . OD 0 
R = 0. 0 DO 
DO 150 J = K,L 
IF (J -EC- I) 



150 



C =C + D A B S ( A ( J , I ) ) 
fi= R + DA ES 
CONTINUE 



GO TO 150 



c 




-GU ASL 




IF (C .SC. 
G= R/ 3 A DIX 


0 . GDO 




F= 1 .00 0 
s=c*a 




160 


IF (C -GE. 
F= F* RADIX 


G) GO 




C =C *92 
GC TO 160 




170 


G =R * 3 A DI X 




180 


IF (C .LT. 


G) GO 



) C OR 
0.0 DO) 



R DUE TO UNDERFLOW-- 
GC TO 220 



GO TC 170 



190 



109 



c 

190 



200 

210 

220 

230 



F»F/aADIX 
C=C/32 
GC TO 180 



IF ( (C + 

G-1 .0D 0/F 



B) / F . GE. 



NOW BALANCE — 

0.9 5 CO * S) GO TO 



220 



SCALil(I) =SCALE (I) *F 
NCCON V ~ » THUS. 

DC 200 J = K r N 

*$«:£’ •= 

I? (NOCONV) GC TO 140 

L C W = K 

IG H = L 

HETU3N 

END 
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c- 

10 



20 

C- 



30 



40 

50 

C~ 



60 



70 

30 



90 

100 



SUBROUTINE 0RTEE5 ( N 2 , H , LOW , IGH, A ,OR 1) 
INTEGER I,J,2 , N.II , JO, LA , 2P , 21 2, IGH, K PI, LOW 
REAL*B A (fifljjl) # «5RT (IGfi) 

REAL*9 ? ,G,H , SCALE 
REAL* 9 DSQRT, £ AES , DS IGN 
L A - EG it - 1 
KP1=L0Wf 1 



IF {LA .17. K ? 1 ) GO TO 100 

DO 90 2= K PI , L A 

H = 0 . GD 0 

CRT (2) -0 . ODO 

SCALE= C. CDQ 

-------SCALE CCLU2N (ALGOL TO L THEN NOT NEEDED) 

DO 10 1 = 55. IGH 

SCALED SCALED CAES ( A ( 1 , 2- 1 ) ) 

IF (SCALE .SQ. C.CDO) GO TO 90 

ai*a+iGH 

DO 20 11=2, IGH 

I=2P-II 



031(1) = A JI,2- 1 ) /SCALE 
H=H*0&? ( I) *0 RT (I) 

CONTINOE 

G = - QS IGN (CSQST (H) ,031(2) ) 
H = H-OR T (2) *G 
05T (2) = 0 FT (2) -G 



DO 50 J = 2,N 
F=0.0D0 

DO 30 II=«,IGH 
I- 2 P- 1 1 

F = F + OR T (I) *A (I ,J) 
CONTINUE 



•FOR 2 



( I - ( 0 * U T ) /H) 



A 



F = ?/H 

DC 40 I=«,IGH 

A (I,J) -A |I,J) - F*CET (I) 

CONTINUE 



FC R2 (I - (U*CJT) / H) *A* (I- (U*UT) /H) 

DC 30 1=1, IGH 
F = C . OD 0 

DC 60 JJ = tt, IGH 
J=2P-J J 

F= F ♦CRT ( J)*A (I,J) 

CONTIN UE 
F=E/H 

DO 70 J=fi,IGH 
A (I,J) = A il, J) - E*0 ET ( J) 

CONTIN US 



GET (2) =SCALS*CST (2) 

A (2,2- 1) =SCAL£*G 

CONTINUE 

RETURN 

END 
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no 



10 

20 



30 



40 



50 

60 

70 

80 



SUE30UTINE ORTHAN (N K, N , LOW , IGH, A , OR T, Z) 

INTEGER I,J, N , KI MP, NH,IG H ,LCW ,MP 1 
3EAL*8 A (SI, I G 5) , CRT (IGH) ,Z (NM,N) 

REALMS G 

C INITIALIZE Z TO IDENTITY MATRIX 

DC 20 1=1, N 
DC 10 J = 1 , N 
Z (I,J) = 0 . CDO 
Z (1,1) =1 .GDO 
C CNTEN UE 
KL=IGH-LCW-1 
IF jKL .IT. 1) 

DC 70 MM=1,KL 
HF = IGII-MM 



GO TO 80 



IF (A<MP,MP-1) . EQ. 0. 0 DO) GO 10 70 

M?1=M?+1 

DC 30 1 = MF 1 , I G H 

ORT (I) =A (I,a?-1) 

DC oO J = MF, IGa 
G =0 . OD 0 

DC 40 I=MF, IGH 
G=G+03 T (I) *2 (I , J) 

DIVISOR BELOW IS NEGATIVE OF H FORMED IN OP.THSS . 

DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW 

G= (G / OFT(MP) ) / A (MF , MP- 1 ) 

DO 50 I=MF,IGH 
Z {I, J) =Z (I, J) ♦ G^O ST ( I) 

CON TIN 02 
CONTINUE 
RETURN 
END 
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SUBROUTINE HQ 3 2 ( N M , N-LOW , I G H , H , W R , W I,Z, ISRH) 

INTEGER I,J, K , L, H , N ,EN, II , JJ , LI, MM, N A, SM, NN, IGH, ITS, LOW, MP 2, SNM2 ,1 
1 EBB 

REALMS H (NM, N) , rf B (N) , VI ( N) , Z ( N S , N) 

3 EAL*8 ?,Q,R,S,T,»,;<,Y,RA,3A,7I,7R,Z2,NCRM,MACHE? 

REALMS DSQItT , CAES, CSIGN 
INTEGER MI NO 
LOGICAL NCTLAS 
COMPLEX *1623 
COMPLEX *16DCMPLX 
REALMS D REAL . DIM AG 

C STATEMENT FUNCTIONS ENABLE EXTRACTION OF REAL AND IMAGINARY — 

C PARTS OF DCUELE PRECISION COMPLEX NUM32RS 

DBEAL ( 23) =23 

DIMAG (23) = {0. 0 D 0 . - 1 . C DO ) * Z3 
DATA M ACHE?/ 234 IGCCOCOOUOCOOO/ 

IERR=0 
NCRM=O.OCO 
K s 1 

C STORE FOOTS ISGIATED 3Y 3ALANC AND COMPUTE MATRIX NORM 

DC 20 1=1 , N 
DC 10 J=K,N 

10 NORfl=NOBM+DABS (H (I,J)) 

IF (I .GE. LO Vi .AND. I . LE. IGH) GO TO 20 
W B (I) = H (1,1) 

VI (I) = G. CDO 
20 CONTINUE 

EN =IGH 
T=0. ODO 

c SEARCH FOR NEXT EIGENVALUES 

30 IF (EN . LT. LCV) GO TO 290 

I TS=0 
N A = EN- 1 
ENM2 = N A- 1 

C LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT 

40 DO 50 LL=LOW,EN 

L= EN+LOW -LL 

IF (L .SC- LOV) GC TO 60 
S = D ABS (H (L-1, L-1) ) +CABS (H (L , L) ) 

IF (S .EC- O.OCC) S= NOR M 

IF (DABS (H(L, L-1) ) . LE. MACHEP * S) GO TO 60 

50 CONTINUE 

C FORM SHIFT 

60 X = H(SN,EN) 

IF (L . EC- EN) GO TO 220 

Y = H (NA , NA) 

W=H (EN , N A) *H ( N A , EN) 

IF (L .EC. NA) GC TO 23 0 
IF ITS . EO. 30) GO TO 500 

IF ITS • NE • 1C .AND. ITS .NS. 20) GC TC SO 

C FCR tt EXCEPTIONAL SHIFT 

T = T+X 

DC 70 I=LCW,ZN 
70 H (I f I) «H (1,1) -X 

S = CABS (H (EN, N A) ) ♦EAESfH (NA, ENM2) ) 

X =0. 75 DO * S 

Y = X 

V=-0.4375D0*S*S 
30 ITS=ITS + 1 

C LOCK FOR T VC CONSECUTIVE SMALL SUB-DIAGONAL ELEMENTS. 

DO 90 3 M = L , EN M 2 
H=ENM2+L-MM 
22=H (M , M) 

R=X-ZZ 
S = Y-22 

P*(H * S - V) / H (M+1 , M) ♦ H ( H r M ♦ 1 ) 

Q=H ( M+ 1, 3+1) -2Z-R-S 
R = H (M+2 , a+1) 

S = D ABS (P) + DA3S (Q) +CAES (R) 

P*P/S 

Q=C/S 

a=a/s 

IF (M .EC- L) GO TO 100 

IF (DABS (H(M, M-1) ) * (DABS (Q) + DABS (8)) .LE. MACHEP * DABS(P) 
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90 

ICO 



1 10 
C— - 



120 



130 

140 



150 
160 
C 



170 

180 
C— - 



190 

200 

210 



C 

220 



C — - 
230 



1 * (DA3S (H(M-1,H-1) ) + 0 A 25 (2 2) + DA ES ( H ( 3 + 1 , a H ) ) ) ) GO TO IOC 

CONTINUE 
«P2=1+2 

DC 110 I=ttP2,£N 

H (I , I- 2) =C. OD 0 

IE (I . iC- 3P2) GC TO 1 10 

a (1,1- 3) -C.OD 0 

CONTINUE 

DOU3LE QR STEP INVOLVING ROWS L TC EN AN 0 COLUNNS a TO EN — 

DO 210 K s«,aa 
NCTLAS-K.SE. A A 
I? (K .EC. H) GO TO 120 
P*H (K, K- 1) 

C=H(K* 1,K-1) 

8 = C . OD 0 

IE {NOT LAS) 3= H { K +2 , K- 1) 

X=DA3S (P) + QA3S JQ) + C AES (3) 

IE (X . EC- O.OCO) GO TO 210 
P = P/X 
Q=C/X 
R ®R/X 

S = CSIG N { ESQRT (P*F+Q*C+H*B),P) 



IF (X .EC* 2S) 
H (K,X- 1) = -S*X 
GO TO 14 C 



GO TO 130 



IE {L 

P-P+S 

X=F/S 

Y*Q/S 

Z Z = 8/S 

Q=Q/P 

R= S/P 



. X E • a) B (K,X- 1) =- H (K, X- 1) 



■SOW aOCIFICATICS- 

DC 160 J = K,N 
P = H (K, J) + C*H (K + 1 ,J) 

IE (.NOT. NOTLAS) GO TO 150 



P=P+a*H ( K + 2 . J) 

H (K + 2, J) =H(K+^,J) -?*ZZ 
H (K+1 , J) *H K+'l ,J) -P* Y 



H ( K , J) * H <K, J) -F*X 
CONTINUE 
J -flINO ( EN ,K +3 ) 

COLUMN aCCIFIC A HON- 
DO 180 1*1, J 
P-X* H ( I , K) + Y * B (I , K + 1 ) 

I? (.NOT. NOTLAS) GO TO 170 
P=P+ZZ*H { I, K + 2) 



H (I, X+2) *K (I , K + 2) -P*R 
H (I , K + 1 i = H(I,K*1) -P*Q 
Hjl ,X) *H (I , K) -P 
CONTINUE 



, (I, K+2) =Z(I , K + 2) -F*B 
! (I,K+ 1) =2{I r K + 1 ) -P *Q 



ACCUaULATE TRANSFORMATION S- 

DO 200 I = LOW r IGH 
P=X*Z ( I, K) + Y*Z (I, K+1) 

IF (.NOT. NOTLAS) GO TO 190 
P=P+ZZ *Z {I, K + 2) 

zr ~ 

z 

Zjl.Kl -Z (I ; K) -P 
CONTINUE 
CONTINUE 
GO TO 40 

0'IE SCOT FOUND— 

H (EN. EN) = X+ T 
W3 (EN) = H { EN , EN) 

HI (EN) =0 . CDO 
EN-N A 
GO TO 30 

T»C ROOTS FOUND — 

P= (Y - I) /2 . OD C 

C=p*p+u 

ZZ-D SQ R T (DABS (C) ) 

H (EN , E N) = X +T 
X = H (EN ,EN) 

H (N A, N A) - Y+T 

IF (Q .LT. 0.0D0) GO TO 270 
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■SEAL PAIR- 



240 
C 



25 0 
C~- 



260 

C 

27 0 



280 



C — 
C — - 
290 



C 

30 0 



310 

320 



330 



C 

340 



HB (EN) =X -ff / 2Z 



zz = ?+dsign<zz,p) 

HR (N A) = X + ZZ 
HR (EN) = HR (N A) 

I? (ZZ . NE . 0.000) 

HI (NA) = 0 . 0 DO 
HI (EN) = C . CDO 
X= H (EN , N A) 

S = D AS S (X) +DA3S (ZZ) 

P=X/S 
Q = Z Z / S 

R=CSQRT (P*?+Q*C) 

P = P/R 
C=C/H 

ROW flOD IF ICATI CM 

DC 240 J=NA,N 
ZZ=H ( N A , J ) 

H (NA,J)=C*ZZ + P*R(EN,J) 

COLUMN MCDIFIC ATICN 

DO 250 1=1, EN 
ZZ= H (I , N A) 

H ( I , N A ) =C*ZZ*F*H (I, BN) 

H ( I , EN ) = C *H ( I , £N) -P*ZZ 
CCNTIN UE 

ACCUMULATE TRANSFORMATIONS 

DC 260 I=LOW, IGE 
ZZ = Z (I , NA) 

Z ( I , N A ) =C*ZZ + F*Z (I, EE) 

Z(I,ZN)=C*Z(I,EN) -F* ZZ 

CONTINUE 

GC TO 280 

COMPLEX £ AIR 

HE (NA) = X + P 
HR (EN ) = X + F 
HI (NA) =ZZ 
HI (EN) =-22 
E N= ENM 2 
GC TO 30 

ALL FOOTS FOUND. z ACKSU0 S TIT UTE TO FIND- 

VECTORS C? UPPER TRIANGULAR FORM 

IF (NORM .ZQ. 0.010) GO TO 510 
DO 450 N N = 1 , N i 

EN=N+1 -N N 
P = HR (EN) 

Q=WI (EN) 

N A =EN- 1 

IF (Q) 370,300,450 

SEAL VECTOR 

M= EN 

H (EN , E N) = 1 . QD 0 

IF (NA . EC. 3) GO TO 450 

DO 360 11=1 , NA 

1= EN- 1 1 

H = H (1,1) -P 

a=H (I, EN) 

IF (M .GT. NA) GO TO 320 

DO 310 J = M , N A 

R=R + 9 (I, J) *9 (J ,2N) 

IF (Hi (I) .G*. 0.0E0) GO TO 330 

ZZ = H 
S = B 

GO TO 360 
«=I 

IF ( HI (I ) .NE. 0 . CDO ) GO TO 340 
T = W 

IF (H .EQ. 0.0D0) T=MACHEP* NO Eff 
H (I,2N) =-B/T 
GO TO 360 

SOLVE REAL EQUATIONS 

X = H (I, 1+ 1) 

T = 9 (1+ 1,1) 

Q = (HR ( I ) - P )* (HR (I ) - ?) ♦■HI (I) * WI (I ) 

“ - ZZ * R ‘ 



T= (X * S 
H ( 1 , 2 N ) = T 



R)/C 
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350 
360 
C — 



c — 
370 
C — - 

c— 



IF (DABS (X) 

a<i*i, 23 n«< 

GC TO 360 
H(IM,EN) = (-S - Y 
CONTINUE 

GO TO 450 



.IE. C A3 S (ZZ) ) 

a - a * i)/x 



GO TO 350 



i)/zz 

END SEAL VECTO 5 

CO 3PL EX VECTOR 



380 

390 



400 



41 Q 



C 

4 20 



430 



440 

C 

45 0 
C 



460 
470 
C — • 
C— 



LA SI V -CTC3 C0.1PONF5T CHOSEM 13 AGIN AST SO 

EIGE3V SCTG3 J3AI3IX IS TP. I A N G U L A ct 

IF (DABS (H (Ell , SA) ) . LE. DABS (H |NA,2N ) ) ) GOTO 380 

H (NA, N A) =Q/H ( E N , S A) 

G<! :, r6 3 3?3<r (a (EN ' £N) * ?)/a(EH,NA) 

23 = DC.1 PL X |0.3D0,-H(NA,EN) )/DC IELX (H ( NA, SA) 2) 

H (HA, HA) =C32AL (Z3) 
a (NA,2N) =DIHAG (23) 
a (E'l , H A) =C.0DC 
H £N, EN) =1.0 DO 
2Na2=N A- 1 

IF (EH. 12 .2Q. 0) GO TO 450 

DO 440 11=1, ENE2 
I=N A-I I 

S A =H (I r EN) 

do 400 j=«,:ia 

3 A = 3 A+ H (I, J) * H ( J r NA) 

SA=SA+H (I , J) (J, ENJ 
CONTINUE 

IF ( W I ( I ) .GE. 0.0C0) GO TO 41C 

zz=w 

5= B A 

S = 

GC **?0 44 0 
M = I 

IF (Wifi) . N E . C.OCO) GO TO 420 
Z3=DCMPLX (-3A ,-SA)/DCMPLX (W , Q ) 
a ( I , N A ) = CHEAT (Z3) 

H il.SN) = C I il A G (23) 

GO TO 440 

SOL7S COMPLEX EQUATIONS 

X=H (I, 1+ 1) 

Y=H I* 1 ,1) 

7 R= (WR (I) - ?)*(WP(I) * P)+WI (I)*WI(I)-Q*Q 
71= (VB (I) - ?) *2 . 0D0 *Q 

I? (VH .FQ. O-OLJ .AND. 71 .SC. 0.0D0) 

1 A BS (Q) ♦ LABS ( X) * CABS(Y) ♦ EABS(ZZ)) 

Z3=DCM PL I (X*R- ZZ*RA+Q*S A , X*S-ZZ*SA-Q*RA)/DCMPLX (VR,7I) 
H (I, NA) = CHE AL (Z3) 

H JI, EN ) = CIMAG (23 

IF (DABS ( X) .IE. CABS (ZZ) ♦ DAES(Q)) GO TO 4 30 
H (1+ 1 , N A) = (- 3 A - * * H ( I , N A ) ♦ Q * H(I,EN))/X 
H (1 + 1 r EN i = (-S A - i * H I,SN) - Q * H(I,NA))/X 
GC TO 440 

Z3 =DCM PL X (-2- Y*8 (I , N A) , -S -Y * H |I,SN) ) / DC H P LX ( ZZ , Q) 

H (1+1, NA = DR E AL (Z3| 



:hat- 



VR= MA CH EP*NO RM* (DABS (*) > D 



H (1*1, EN ) -DI tt AG ( Z 
CONTINUE 



END COMPLEX VECTOR — 

SUBSTITUTION. 7ECTORS OF 



ISOLATED HO DIS- 



CONTINUE 

END EACK 

DO 470 1=1, N 
IF (I .GE. LOS .AND. I . LE. IGB) GO TO 470 
DO 460 J = 1 , N 
Z <I,J) - H (I, J) 

CONTINUE 

MULTIPLY EY TRANSFORMATION MATRIX TO GIVE- 

7ECTCES OF ORIGINAL FULL MATBIX. 

DO 490 J J=LOW , N 
J = N + LO Vi- J J 
a=aiN0 (j , igh) 

DO 490 I =LCS , IGH 
ZZ = 0. 0 DO 
DO 480 K=LOW,3 
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480 ZZ =ZZ* Z (I ,K) * H { ? f o) 

Z (I, J) = ZZ 
490 CONTINUE 
GC TO 510 

C SET ERROR ~>NC CONVERGENCE TO AN 

C EIGENVALUE AFTER 30 ITERATIONS 

50 C IE?R=SN 
510 E ETC RN 
END 
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UUu 



c== 



10 

20 

30 



40 

50 

60 



SUBROUTINE 3 ALEAK ( N K, N f LOW f IGE f SCALE, 

I NT E 32 a I,J,K,D,N,II r NM,ISH,LCW 

REAL** SCALE (N) ,Z , S 

I? . EQ . 0) GG 1C 60 

IF (IGH .IQ. LOW) GO TO 30 

DO 20 I=LCW,IGH 

S=S CAL E (I) 



LEFT HAND EIGENVECTORS ARE BACK TRANS FOR a ED 

1 ? x 5 E FOREGOING STATEMENT IS REPLACED 3 Y — 

S^r 1 „ CD C/ SCALE (I) . 

DC 10 J=1,« 

uh i\n a ' j) * s 

do so ii=i, y 
i=i i 

I? (I -GE • LOS -AND. I . LE. IGS) GO TO 50 
IF (I . LT. LOW) IsLCW-II 
X=SCAL2 (I) 

IF (K A C- I) GO TO 50 
DC 40 j=i,a 
S = 2 (I, J) 



Z (I, J) =2 ( K , J) 
2 ( K , J ) =S 



CONTINUE 

CONTINUE 

RETURN 

END 
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SUBROUTINE HQ R (NE , N ,LOW , IGH , ?.,3R ,31 , IEFR) 

INTEGER I f J, K , I f .%3 # ZN, LL ,A M , N A , NM , I GH, ITS , L03 , A92 , 3NH 2 , IE HR 
HEAL*8 H (NM,N) ,*5 (N) ,WT ( N ) 

REALMS P,C,fi,S,T,W ,X,Y,ZZ,NORE,MACHEP 
REALMS D SORT, CAES ,CSIGN 
INTEGER KINO 
LOGICAL NCTLAS 

DATA MACE EP/Z 3 il 10C0CCQQQ0QQ00/ 

X =0 

ncb8=o. oto 

Xsl 

STORE FOOTS ISCLAIED BY BAIANC AND COMPUTE MATRIX NORM 

DO 20 1= 1 ,3 
DC 10 J= K N 

NORM=NORM*DA3S (H (I,J)) 

I? (I .GE. LOS .AND. I . LZ. IGH) GO TO 20 
WH(I)*H(I,I) 

HI (I) =0. CDO 



C — 
10 

20 



c — 
30 



C — 
4 0 



50 

C 



70 



30 
C — 



CONTINUE 
E3«IGH 
T = 0. 0 D 0 



-SEARCH ?0 R NEXT EIGENVALUES 

IP (2N .LI. LOS) GO TO 250 

ITS=0 

NA=ZN- 1 

SNM2=NA-1 

.11 1 LOCK PCS SINGLE S M ALL SUB-DIAGONAL ELEMENT- 

DO 50 LL = LOW , Z N 
L= EN+LON-LL 

IP (L .EC. LOS) GC TC 60 
S - CABS (H (L-1,L-1) ) ♦DABS (H (L,L) ) 

IF (S .EC. O.OCG) S= NOR M 

IP (DABS (H(L,L-1) ) . LE. 3ACHEP * S) GO TO 60 

CONTIN CE 

•FORM SHIFT 



0 X = H | 


[EN , EN) 


I? < 


L . EC. 2N) GO TO 


Y = H { 


[N A ,3 A) 


3=H | 


[EM, 3A)*H ( N A , E N) 


IF { 


;L .EC. Nrt) GC TO 


IF ( 


ITS . EO. 30) GO 


' IF { 


llTS • NE • 1C .AND 



90 

100 



110 



FORM EXCEPTIONAL SHIFT 

T = T ♦ X 

DC 70 1= LC 3 , ZN 
H (1,1) =H (1,1) -X 

S = CA3S (H (EN, N A) ) ♦ CASS(H(NA,ENK2) ) 

X=C.75D0*S 

Y=X 

W=-0.4375D0*S*S 
ITS=ITS+ 1 

LOOK FOR TWO CONSECUTIVE SMALL SUB-DIAGONAL ELEMENTS. 

DC 90 MM=L,2NM2 
«=ENM2 +L-MM 
ZZ = H (M , M ) 

R=X-ZZ 
S= Y-ZZ 

P = (R * S - W) /H (M + 1 ,M)+H (E,M*1) 

Q=H (M + 1, M + 1) -ZZ-R-S 
R=H (E+2, M+1 

S = D A3S (P) ABS (C) ♦CAES(S) 

P = P/S 
Q=C/S 
R = R/S 

IF (M .EQ. L) GO TO 100 

IF (DABS (H M , 3 - 1 ) ) * (DABS (Q) ♦ DABS (R) ) .LE. 3ACHEP * DABSf?) 

* (DABS (H (M- 1 ,M- 1) ) ♦ DAES (ZZ) + DA ES ( H (M «■ 1 , M ♦ 1 ) ) ) ) GO TO 100 
CONTINUE - 

M F 2= M ♦ 2 

DO 110 I =MP2 # F N 

H (1,1-2) =0.0D0 

IF (I .EQ. 3 P 2 ) GC TC 1 10 

H (1,1- 3) =0. 0 D 0 

CONTINUE 
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120 



130 

140 



150 

160 

C 



170 

130 

190 



C — - 
200 



C 

210 



•-D0U3LE QR STEF INVOLVING SOWS L TO EN AND COLUMNS 
DC 190 K=M,NA 
NOTLAS = K . NE . M A 
If (K . SC • M) GO TO 120 
p=y (K, k- i) 

Q=a (K+ 1, K- 1 ) 

H =C. 0 Q 0 

IF (NOTLAS) R = E ( K ♦ 2 , K-1) 



TO 



I -DA3 3 (P) +0A3S (C) *D AES (H) 
0 . 0 CO) “ ““ - 



GO TO 1 90 



IF (X . EC- 
P=P/X 

c=c/x 

a=a/x 

S = CS 13 N ( CSQRT ( P * P *C *C+R * R ) , ?) 



IF (K .EC. M) 
H (K, 5- 1) = -3*X 
GO TO 14 C 



GO 10 130 



[f (L 
P = P+S 
X = P/S 
I=Q/S 

zz=r/s 

Q=Q/P 

a=a/? 



N E . M) H(K,K- 1 )»-H(K,EC- 1 ) 



DO 160 J=K, ZN 



•ROW MODIFICATION- 



P=fl (X, J) + Q *H ( K * 1 ,J) 

IF (.NOT. NOTLAS) GO TO 150 
P=P>R*H (K*2, J) 

3 (K + 2, J) =H(K+2,J) -?*ZZ 
H (X+ 1 , J) -fi ( K ♦ 1 ,J) -?* ¥ 

H (K,J) = H (K, J) -F*X 

CONTINUE 

J=ttIN0 (EN,K*3) 

COLUMN MODIFICATION- 

DO 180 I =L , J 

P=X*H (I,K) + Y * H ( I , K ♦ 1 ) 

IF (.NOT.* NOTLAS) GG TO 170 
P=P+ZZ*H (X,K + 2) 

H (I,K*2) =H (I , K+2) -P*R 

H (I,X + 1) =8(1, K ♦ 1 ) -P*Q 

H <I,X) -H (I , K) -E 

CONTINUE 

CONTINUE 

GO TO 40 

ONE ROOT FOUND 



C 

220 



230 

C 

C 

240 

250 



UR (EN) = X ♦ T 
UI (EN) =0. ODO 
EN^NA 
GO TO 30 

TWC hqcts FOUND- 
RY (Y - X) /2- 0C0 

Q-p*p+ 

ZZ=DSQHT (DABS (C) ) 

X = X+ T 

IF (Q .LT. O.OCO) GO TO 220 

REAL f MR 

ZZ = P+DSIGN(ZZ,P) 

UR (NA) =X + ZZ 
STB (EN) = W E (NA) 

IF (ZZ .NE. O.ODO) UR (EN) =X-W/ZZ 
UI (NA) =0 . ODO 
W I (EN) = 0 . ODO 
GC TO 230 



UH (NA) =X ♦F 
UR (EN) = X + P 
HI (NA) =ZZ 
HI (EN) =-ZZ 
EN= ENM 2 
GO TO 30 



-COMPLEX P AIR- 



IER S = E N 

RETURN 

END 



-SET ERROR — NO CONVERGENCE TO AN- 
LGENVALUE AFTER 30 ITERATIONS 
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nnnnonnnnonnnnn 



SU3H0UTINE PSDCAI (N2 ,NS . FA , X ^C , GS , G7 , C , UC , HT, HO , H , 

1 FBGE , NG,GAM , ACL , i , WR, 'i 1 ,0 1 , CZ,JCF, BBS , i, 3 . 3B ,CC , X XU , 

2 I PSD , INOEM ) 

« PSDCAI COMFCIES IKE PSE CP OUTPUTS CE CONTROLS OF 



CONTROLLED 


SYSTEX 






= 


IYU= 1 


CUT PUT PSD 






= 


= 2 


CONTROL PSD 






— 


= 3 


EOT H OUTPUT AND CONTROL PSD 






2 


IPSC=1 


FSD 






= 


=2 


PSD AND TF RESIDUES 






- 


I N C B X = 


1,2,.-. NG NORMALIZED 3Y I Iff 
KG + 1 , • • • NG + l> C NORMALIZED BY 


PROCESS 


NOISE 


I 




ITH MEAS 


NOISE 


— 



DOUBLE PRECISION FA. X ,G W .GV ,C EG E , G AM , ACL, F, W S , 31, D 1 , £2 , RES, 

1 EE.CC, Q, B.PSC, W, CNORt!,DN 1 ,E M A X ,ZLOG , EMOD , D h , ST ,0M , A £ , A I , H U , CW 1 
COMPLEX * 1 6 ZD 2N ZZ 

DISEHS ION PA (N2 A , N2) , X (N 2, N2) . GW (N2,N G) , C (NC,NS) ,H Y ( NO , N2) ,H ( NC. NS) 

1 , F3GI ( NS , NO) ,G AX (NS , NG) , ACL( NS, NS) , F ( NS , NS) , W B ( N2) , WI(N2) , Dl (N2) ,D 
22 ( H2) , 3£S (N2) , C (N G. N G) ,3 (NO , NC) # PSD ( 30) , W (30 ) , 3S (:i 2 ) , C C (N 2 ) , G V ( N 2 , 
3 NC) , Htl ( NC , N2) , EH 1 (4) 

INTEGER JC? ( N 2 ) 

DATA DH1/1.D0,2.DC,5.D0, 10. DO/ 

I? (TYC . EQ. C) I Y U= 1 
IF (I^OSM • EQ • 0) INORil=1 
IPT=0 

IF (I? S D .GT. 1) IPT=1 
IX=IN0 3K-t!G 

IF (IX . GT. 0) WRITE (6,330) IX 
I? (IX .IE. 0) WRITE (0,340) INCRX 
N S Q= N 2 * N 2 

CD XP GTE EIGEN SYS TEN OF CONTROLLED SYS TEX ; FORM FA 

CC 10 1=1, NS 

.DO 10 J= 1 ,N S 
FA (I, J) =ACL (I , J) 

FA (NS+I,J)=0.D0 

DO 30 1=1, NS 

DO 30 J = 1 , N S 

ST=0. DO 

DO 20 K = 1 ,NO 

ST=ST+ FBGr (I , X) *H (K, J) 

FA (I, MS* J)=-ST 

FA (NS + I, HS+J) =F (I, J) -ST 

CALL RAP ENT ( N 2 , N 2 , M 2, 9 , F A, 4 , ■ ( 9 ( IX , 1PD 1 3 . 6) ) • ) 

DEBUG ABOVE 

CALL 3ALANC ( N 2 , N 2 , F A ,LOW , I HI GH, D 1) 



10 

20 

30 

C — - 



CALL OBTHES ( N2 , N2 ,ICW, I HIGH , FA, C2) 
CALL ORTF.AN l N 2, N 2 , LOW, I B IG H , f A , D2, X ) 
CALL HQR2 (N 2 , N2 , LO W ,IH I G H , F A , W R , W I , X , 
IF (IERB .NE. 0) GO TO 320 



IFSR) 



CALL 

CALL 



40 

50 



C — 
60 



70 

80 



3ALBAK (N2,N2 , ICS. I HIGH. C1,N2,X) 

RAP ENT (N2,N2,N2,9, X,4, r (9 (IX, 1PD13.6)) *) 

DEEUG ABC YE ; DETERMINE XCDAL X AT BICES 

IF (ITU .EQ. 1) GC TO 60 

2 i HSUBU 

DO 50 1=1, NC 

DC 50 J = 1 , N 2 

ST = 0 - 0 0 

DC 40 K= 1 , NS 

ST = ST-C (I ,K) *X <F ,0) 

HO (I, J) = ST 
GO TO 90 

HSUBY 

DO 30 1=1, NO 
DC 3Q J = 1 , N 2 
ST =0 . D 0 
DC 70 K= 1 ,NS 

ST=ST+H (I ,K) * X (F,0)-H (I , K ) * X ( KS + K , J) 

HY (I, J) = ST 

CALL HAP RKT (NC,NC,N2,9,HY. 4, * (9 (IX, 1PD13.6) ) ') 
DEBUG AEOVE 
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on 



90 



100 
1 10 



120 

c— 



1 ) 



130 

140 
C — 



150 



160 
C— - 



170 
C — - 



CALL MINV (NSQ,X,N2,ST, D1 ,02) 

CALL a&PHST (N2,N2,N2,9 ,X,4, • (9(1X r 1PDl3.o)j ») 

— — — — —————— 0 E DO G A 50 V 2— • ——— — —————— — — —— — - 

GSUB W 

DO 110 1=1, M2 

DC 110 J=1,NG 

ST = 0 . 0 DO 

DC 100 K = 1 ,MS 

ST=ST-X (I, NS *K) ♦GAM (K,J) 

g s (I , J) -ST 

CALL a APS NT (M2,N2,MG,9,GW,4, 1 (9 (IX, 1 ? D 1 3 . 6) ) # ) 

DEBUG ABOVE; USE SELECTED NORMALIZATION 

IF (IN OR a .LS. NG) DNORM* 1. DO/C (INOR H,I NOSH) 

IF (INORM .GT. NG) DNORM= 1. DO/R (INCH M-iiG, INORM-NG) 

DETERMINE BANDWIDTH Gr CONTROLLED SYSTEM- 

EM AX =0 . DO 
DO 120 1=1, N2 

EMCD=0 AES (S3 (I) **2 ♦ WI(I)**2) 

I? (EMOD •GT. EMAX) EMA X = EMOD 

CONTINUE 

2 MOD=D S Q ET (EM AX) 

EMO D=2 *E MCD 

SOUND OP TO NEAREST 2,4, 5,6,10 

ELOG = D LOG 1 0 ( E MOD ) 

IF (ELOG . Lt . O.CC) IPO W = -ID I NT (DAES (ELOG) 

I? jELOG .G2. O.DC) IPO S= ID I NT (ELOG) 

2MAX=EMOC*10** J-IPOfc) 

I? (EMAX .GT. 2. EC) *KOD = 2.D0 
IF (UMAX .GT. 4. DC) ZM0D=4. DO 
IF ( EM AX .GT. 5. DC) EMOD = 5.DO 
IF (EMAX .GT. 6. DC) EMOD = 8.DO 
IF (EMAX .G2. 10. TO) EMOD = 10.D0 

S M A X= E M 0 D * 1 0 * * I E C * 

DW= EM AX/20. DO 

ADD POINTS 3 DECADES UP- 

IF (EMOD .LT. 5.0) GC TO 130 
E M A X = 1 . OD 1 
IK=3 

GC TO 1 4 C 
EM AX=5 .DO 
IK =2 

CONTINUE 

STCRE 3 0 FR ECUENCI ES 

DC 150 1=1,20 
a (i) -os* (i-i) 

DC 160 1=1,3 
IP=20* 3 * (1-1 ) 

DO 160 J=1,3 
IX=MOD (IR + J-1 ,3) +1 
J J-0 

IF (IK .EQ. 2 .AND. J .GE. 2) JJ= 1 
H (IP+J) = CW1 (IX) *10** (IPOW+I- 1 +JJ+IK- 2) 

CONTIN UE 
IX = MOD (IK, 3) + 1 

W (30) =D W1 (IX) *10** (IPOW+ 3 +IK-2) 

LARGE LOOP THRU OUTPUTS — - 

IF (IYU .ZQ. 1) NI=NC 
IF (IYU .EQ. 2) NL= N C 
DO 310 L = 1 , N L 
DO 170 1=1,30 
PSD (I) =0 . DO 



DO 220 1*1, NG 
D N 1= ON CRM*Q ( I , I) 

1 “ “ 



•LOOP THRU PROCESS NOISE- 



(IYU 


• z, Q . 


1 .AND. IPT .EQ. 


1) WRITE 


(6,350) 


(IYU 


• LQ. 


2 .AND. IPT .EQ. 


l[ WRITE 


(6,360) 
, G * , N L , 


ITU. 


HQ. 1) 


CALI RESID (I,L,: 


32, JC? , NG 



180 



IRES, 33.CC ,I?T) 

IF (IYU. EQ. 2) CALL RSSID 
IRES. 33 , CC .IPT) 

DO 210 K= 1 , 20 
ZZ=DCMPL X (0 . D 0 , 0 . DO) 

OM=W (K) 

DO 200 11=1, N2 

IF (WI (II) ) 20C , 180, 190 

ZD=DCMPL X {-SR (II) ,Cfi-«I (II) ) 



(I, L ,N2 r JCF ,NG ,GW,NL, HU , SR ,SI , 
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190 



200 
210 
220 
C 



230 

240 



250 

26 C 

27 0 



280 

290 

300 



310 

320 



ZZ=RES (I I } /ZD ♦ ZZ 
GC TO 20 0 
RF = *R (II) 

AI=WI (II) 

ZD=DCM?L2 (RB**2 * AI **2 - OM**2,-2. D0*R£*OM) 

zn=dcmflx (res (ii+i) *ai-res (ii) *re,res (i i > *om) 

ZZ=ZZ+ZN/ZD 

CONTINUE 

PSD (X) = ?SD (K) + DN1 * (ZZ*DCONJG (ZZ) ) 

CONTINUE 



— G 303 V 

DO 240 1=1, M2 
DO 240 J = 1 , N 0 
ST=0 . DO 
DO 23 0 K = 1 - N S 

ST=ST*:< ( I ,X) * F EG Z (K , J) + X (I,NS+K) *FBG E (K,J) 

G V (I , J ) = ST 

CALL RAPE NT ( N 2 , N2, NO, 9 f G V, 4 , » (9 ( IX , 1?D 13 . 6) ) * ) 

DEBUG ABC7S, LOO? THRU 3 E AS NOISE 

DO 300 1=1, NO 
DN1=DN0RM*R (I .1) 

IF (IYU .HQ. 1 .ANC. I?T .EQ. 1) WRITE (6,370) I,L 

IF (IYU .Eg. 2 .AND. IPT .EQ. if WRITE (6,380 I,L 

IF (IYU. EC. 1) CALL RESID (I , L , N2 , JC? , NC ,G V , NL , 3Y , A fi , W I , RES , 

1 3B , CC , IPT) 

IF (IYU. EQ. 2) CALL RESID (I , L , N2 , JCF , NO , G V , N L , HU , W R , W I , ? 2S , 

1 3B , CC , IPT ) 

DO 290 K=1.30 

zz=dcm?lx (O.DC,o.ro) 

OK*W (K) 

DC 270 11=1, N2 

IF (WI (II) ) 27 0,250. 260 

Z D= DCM PL X (-WR (II) ,0S-WI (II) ) 

ZZ=ZZ-r RES (II) / ZD 
GC TO 270 
BE = WR (II) 

A 1= WI ( II) 

ZD = DCMPLX (RE**2 «• AI**2 -CM**2 ,-2. DO *RE *OM) 

ZN=DCMPLX (RES ( 11+ 1) * A I- RES (I I ) *RE ,52 S (II ) *OM) 

ZZ=ZZ* ZN/ZD 
CONTINUE 

IF (IYU .EQ. 2 .CH. I .ME. L) GO TO 230 
PSD (K) = PSC(X) ♦ CN 1 

PSDjK) =PSD (K) + D N 1 * (ZZ*DCONJG (ZZ) ) 

CONTINUE 
CONTIN UE 
I? (ITU .EQ. 



) WRITE i 


(6, 390) 


L 


) WRITE 


6 ,400 


I 


(W(I) #SSI 


5(11,1=1,30) 



C 

330 
34 0 

3 50 

360 

370 

380 

390 

400 

410 



CONTINUE 
RETURN 
CONTIN UE 

CALL EREXIT ( N 2, E A , IEHR) 
RETURN 



FOHKAT 

FORMAT 

He 



( / ,4 1 H SUBSECUENT PSD IS NORMALIZED 3Y ME AS NO., 13) 

(/ ,50 H SUBSEQUENT PSD IS NORMALIZED BY PROCESS NOISE NC.,I3 



■’CBMAT ( /38 H TRANSFER FUNCTION FRCM PROCESS NOISE ,I2,3H T0,13H ME 
1 ASU REM ENT ,12) 

FORMAT (/38H TRANSFER FUNCTION FRCM PROCESS NOISE , 12, 3H TO,9H CCN 
1 T RC L 12 ) 

FORMAT (/ 36H TRANSFER FUNCTION FBCtl ME ASU BEM ENT ,12, 16a TO HEASU3S 
1HENT ,12) 

FORMAT (/36H T8ANSFE8 FUNCTION FROM MEASUREMENT ,I2,12HTO CONTROL 

1 , 12 ) 

FORMAT (/1UH PSD CF CUT? UT, 1 3 , 32H FORCED BY ALL NOISE-(3AD FREQ,, 
1 1 5 HNOS M ALIZ2D ESC)/) 

FORMAT ( / 1 5 1! PSD CF CCNT B CL , I 2 , 32H FORCED BY ALL N0ISE-(3AD FBEC, 
1 , 15HNORMALI2EE PSD)/) 

FORMAT (4 (IX, 1H(, El 1.4, IB,, 31 1.4, 1H) )) 

END 
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non 



S 0 ERCU TINE 2RHXIT ( H , A, IE HR) 

EREXIT RETURNS THE NUMBER OF THE EIGENVALUE WHERE HQR2 
FAILS l TH£N_STCP3 THE PROGR AM._ 

I NTEG E5~ IERR 
DOUBLE PRECISION A 
DIMENSION A ( N , N ) 

WRITE (5 1 10) IIP?. 

CALL R A P fa NT ( N , M , N , 9 f A , 4 , 1 (9{ 1 X f 1 ? D 1 3.6)) *) 

RETURN 

10 FORMAT (35H FAILURE IN HQR2 ON EIGENVALUE NO. ,13) 

END 
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c== 

c 

c== 



10 

20 

C~ 

30 

40 



50 



60 

70 



80 



90 

100 

110 



C 

120 

130 

140 

150 

160 

170 

180 



SUEF0U7INE REACF (NS,ISAF,3A) 

INTER ACTIVELY ENTERS _THE "F^ MATRIX ELEMENT BY ELEMENT. 

RE AL*8 EA(NS,NS) .DU M .A NS H 
INTEGER I,J,X,L, i1nS,I3AF 

' ” ■ "Z/'sv 

GO TO 40 



WRITE (5,120) I,J 
CALL BOREAL ( ANS R) 
BA (I,J) = A NSR 
CONTINUE 
CONTIN OE 



data n/'r/ ( i 
I? (IS AF. EQ. 1) 
WRITS (5.130) 
DO 20 i=1, NS 
DC 10 J = 1 , N S 



CALL FRTCMS ('CIRSC5N ') 

CONTINUE 
WRITE (5,140) 

CALL M AT t RT (EA,NS,NS) 

WRITE (5,150) 

CALL R DC EAR (IANS) 

IF ( (IANS. STB. I Y) .AND. (IANS. NE.IZ) ) GC TC 60 
GC TO 70 
WRITE (5,160) 

GC TO 50 
CONTINUE 



IF (IANS. £0.12) 
IF (IANS. SO. IY 
WRITE (5,170) 



GC TC 110 
GC TC 80 



CALL RDINT (IANS) 

K= I ANS 

WRITE (5,180) 

CALL RDINT (IANS) 

L=IANS 

WHITE (5,120) K,L 
CALL BOREAL ( A NSH ) 

DUM= ANS R 
DO 100 I = 1 . N S 
DO 90 J= 1,5s 

I? ( (I .SQ.X) . AND. (J. SQ. L) ) 3A ll, J) = D 0,1 

CONTINUE 

CONTINUE 

GO TO 30 

CONTI SUE 

CALL F ST CMS ( * CL3SCR N ' ) 

RETURN 



FORMAT 
FORMAT 
1SNSION 
FORMAT 
FORMAT 
1 ENT?.// 
FORMA T 
FORMAT 
FORMAT 

1 eIid 



(SX.14HTHE ELEBENT F ( , 12 , 1 H, , 12 , 2H> =) 
(/ ,ix . 36 r. ENT EH THE SYSTEM MATRIX " 



= t STATES NS X * STATES 



X " F" -3 AT EIX ,//, 10X, UIHDia 

. ... .. NS ) 

(// , 1 5 X .33 HTH Z SYSTZ3 MATRIX "F •'-!* AT HI X ...,//) 

{//5x,54aoc yco wish tc change the value of any a at 3 ix elem 

, 10X,t9BTYFE "YES" OR "NO".) 

‘IX, 51 H WARNING: IMPROPER DATA ENTRY! ENTER "YES" OS "NO".) 
5X,5 OH ENTER THE RCW N OSBER OF THE ELEMENT TO BE CHANGED.) 
5X,S3 HENTER THE CCLUHN NUMBER OF THE ELEaENT TO BE CHANGED 
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SUBROUTINE BEACH (NC , NS , 1 5AH , HC) 
INTERACTIVELY ENTERS THE M H" 2 AT El 



IX SEAS ARE HE NT SCALING MATRIX .= 



BEAL*8 EC (NO , NS) ,DUH,ANSR 
INTEGER IANS,I ,J ,8. L ,ISAH 
DATA IY/'Y'/r IZ/ f N*/ 



THIS IS AN EXAMPLE G ? CNE POSSIBLE METHOD CF ARRAY GENERATION 
WITHIN THE PROGRAM USE L?. FOR VERY LARGE DATA ARRAYS, THIS METHOD 
MAY cZ PREFERABLE TO SOME USERS OVER INTERACTIVE ENTRY OF EACH 
INDIVIDUAL ELEMENT. 



Cl 

C2 

C 

C3 
C — 



DO 2 1=1.11 

DO 1 0=1,92 
HC (I , J) 

HC ( 1 , l[ 
HO (2 , 7 5) 
HC (3 , 7 4) 
HC (4,63) 
BC ( 5 , 62 
HO (6 ,76 
HC (7,44] 
HC(3,45 
HC 9, 46] 
HOMO, 47) 
EC (1 1 , 48) 
CONTINUE 
CONTINUE 
GO 10 90 
CCNTINU E 



= O.OD+OO 
= 0.115230 + 00 
= 0.573 CD + 02 
= C. 1 0 00D >0 1 
= 0. 5730D+02 
= 0. 10 00D + 0 1 
= 0.57 3 OD + Q 2 
= C.57 3 CD ♦0 2 
= 0.5730D + 02 
= 0. 57300*02 
= 0.573 OD +02 
= 0.573 CD + 02 



10 

20 

C~ 

30 

40 



50 

60 

70 



IF (IS AH. EC. 1) GO TO 40 
WRITE (5.120) 

DC 20 1=1, NO 
DO 10 J=1 .NS 
WHITE (5,110) I,J 
CALL 3DREAL (ANSR) 

HO (I, J) = ANSR 

CONTINUE 

CONTINUE 



80 

90 

100 



C — - 
110 
120 



CALL F RT CMS (’CLRSCRN ») 

CONTINUE ' 

WHITE (5,130) 

CALL MATPRT (HC,NC,NS) 

WRITE (5,140) 

CALL RDCHAR (IANS) 

IF ( (IANS. N2. IY) .AND. (IANS. NE.IZ) ) GO TO 60 
GO TO 70 
WHITE (5,150) 

GO TO 50 
CONTINUE 

IF (IAKS.ZQ.IZ) GC TC 100 
WHITE (5 , IbO) 

CALL RDINT (IANS) 

K= I ANS 

WRITE (5,170) 

CALL RDINT (IANS) 

L = I A N S 

WHITE (5,110) K , L 
CALL RDREAL (ANSR) 

DU M= AN SR 
DO 90 1 = 1, NO 
DC 80 J = 1 , N S 

IF ((I.EC.X) .AND. (J • EQ. L) ) HO (I , J) = D U M 

CONTINUE 

CONTINUE 

GC TO 30 

COMTIN UE 

CALL F BT CMS ( 1 CLRSCR N ') 

RETURN 



FORMAT (5X,14HIHE ELEMENT H ( , 12 , 1 H , , 12 , 2 H) =) 

FORMAT (/,3X, 5CHENTER THE MEASUREMENT SCALING MATRIX 
1 ,//, 1 0X, 47HDIMENSICN = » OBSERVATIONS NO X 4 STATES 



,, H H - MATRIX . 
NS ) 
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130 FCRNST (//, 13X,4beTHE N5AS0REKENT SCALING NATNIX " H"- NATR I X . . . ,/ 
1 /) 

140 FOSSA! (//5X , 3 4 RCC YCO NISH TC CHANGE THE VALUE OF ANY SATE IX ELIS 
1ENT?,//, 1 OX , 1 9HTYEE "YES” OR "NO".) 

150 F03NAT ( IX, 5 1 HWARNING: ISPSOPER DATA ENTRY! ENTER "YES" OR "NO".) 
160 POASAT (5X,5CHENTES THE SOW NGIBES OF THE ELEN ENT TO 3E CHANGES. ) 

170 F03aAT ( 5X, 5 2 H E ST ER THE CCLU3N NUX3ER 0? THE ELEN ENT TO EE CHANGED 
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nn 



C== 



10 
20 
C — 
30 



40 

50 

60 



70 

80 

90 



S 02 RO UTI N E REACD ( NC , NC , D ) 

ENTERS THE W C W MATRIX MEASUREMENT PEZD-FOEWARD DIST • _ MATE IX 

~iiAL*3~ c7»o7NcrTHu5TANsi”” 

INTEGER I ANS , I ,o . H f L 
DATA II/'Y'/r I2/ f N V 
WRITE (5,110) 

DO 20 1= 1 ,NO 
DC 10 J=1,NC 
WRITS (5,100) I,J 
CALL R DR E AL ( A NS S) 

D (I, J) “A N SR 

CONTINUE 

CONTINUE 

CALL ? RTCMS ( * CLHSCR N *) 

WRITE (5,120) 

CALL M ATF ET ( D ,NC ,NC) 

WRITE (5,130) 

CALL RDCHAR (IANS) 

IF ( (IANS.NE.I X) .AND. (IANS. HE. 12) ) GC TO 50 
GG TO 60 
WRITE (5,140) 

GO TO 40 
CONTIN CE 

IF (IANS.EQ. IZ) GC TC 90 
WRITS (5.150) 

CALL RDINT (IANS) 

K = I A N S 

WRITS (5,16 0) 

CALL RDINT (IANS) 

L= I ANS 

WRITS (5,100) K,L 
CALL RDREAL (ANSI) 

DUM=ANSR 
DC 80 1=1, NO 
DO 70 J = 1 , N C 

IF ( (I. EQ.K) . AND. (J. EQ.L) ) D (I,J) =DUM 

CONTINUE 

CONTINUE 

GO TO 30 

CONTINUE 

CALL FRTCMS ( 1 CLRSCF N ') 

RETURN 



C 

100 
1 10 



120 

130 

140 

150 

160 



FORMAT 

FORMAT 

1 WARD,/, 

2 = #03 
FORMAT 

1 • ,//) 
FORMAT 

IE NT?,// 
FORMAT 
FORMAT 
FORMAT 

1 .) 

END 



i/?5x 

5X , 34 
SER7A 

( // ,5 X , 5CHTH E FEEDFORWARD DISTHIEUTION MATRIX "D” - MATRIX 



4HTEE ELEMENT D ( . 12 , 1H, , 12, 2H) =) 

, 54HSNTER THE MEASUREMENT FEEDTHROUGH MATRIX / FEEDFCR 
H DISTHIEUTION MATRIX "C'-MATRIX . ,//, 8X, 49 U DIME NS ION 
TIOS NC X 4 CONTROLS NC ) 

. £ua 



(//5X 
r 1CX f 
1 X , D 
5X,5 
5X, 5 



, 5 4 B DC YOU WISH TO CHANGE THE VALUE OF ANY MATRIX ELEM 
19HTYFS "YES" OR 11 NC” . ) 

1 ti WARNING: IMPROPER DATA ENTRY! ENTER "YES" OR "NO”.) 
0HENTE3 THE RCW NUMBER OF THE ELEMENT TO BE CHANGED.) 

3 fi ENT HR THE COLUMN NUMBER OF THE ELEMENT TO BE CHANGED 
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SDBROUTI NE READG (NS , NC , ISAG , G) 

ppSACTIVELY ENTSBSJM3 |]G^_2ATRIX_ CONTROL DISTRIBUTION MATRIX = 

SEAL*8 G*<NS,NC) ,DUM,ANSR 
INTEGER IANS, I ,J ,K. L,ISAG 
DATA IY/^'/r 12/' NV 
IF (I SAG . £0 . 1) GO TC 40 
WRITS (5,120) 

DO 20 1=1, !IS 
DO 10 J=1.NC 
WHITS (5,110) I,J 
CALL RDREAL ( A NS R ) 

G (I, J) = A NSH 
10 CONTINUE 

20 CONTINUE 

30 CALL FRTCM5 ( ' CL3SCH N ') 

40 CONTINUE 

WHITS (5,130) 

CALL M A T F 5 T (G ,NS ,NC) 

50 WHITS (5,140) 

CALL SDCHAR (IANS) 

IF ({IANS. NS. II) . AND. (IANS.NE.IZ) ) GO TC 60 
GO TO 70 

60 WBITS (5 , 150) 

GO TO 50 
70 CONTINUE 

IF (IANS.EQ.IZ) GC TC 100 
WHITE (5,160) 

CALL RDINT (IANS) 

K= I AN3 

WRITS (5.170) 

CALL RDlflT (IANS) 

L = I A N S 

WHITS (5,110) K , L 
CALL RDREAL ( A NS 3 ) 

DU AN SR 
DO 90 1=1 ,NS 
DC 30 J=1, NC 

IF f (I .2Q.K) . AND. (J. HQ. L) ) G(I,J)=DUM 
80 CONTINUE 

90 CONTINUE 

GO TO 30 
100 CONTINUE 

CALL F STCMS ( 1 CL3SCHN * ) 

RETURN 

110 FORMAT ( 5 X , 1 4 H TH £ ELEMENT G ( , 12 , 1 H , , 12 , 2H) =) 

120 FORMAT (/ , §X , 5 1 H ENT EH THE CONTROL DISTRIBUTION MATRIX "G" -MATRIX 

1 . 1 0X,43HDIM£NSICN = # STATES NS X * CONTROLS NC ) 

130 FORMAT (//, 1 0 X ,47HTHE CONTROL CISTRIEOTION MATRIX " G" -MAT R IX ..., 

1//) 

140 FORMAT (//5X,54RDC YCO WISH TC CHANGE THE VALUE OF ANY MATRIX SLEM 

1ENT?,//, 1 CX . 1 9HTYFE "YES*' OR "NC".) 

150 FOBMAT ( IX, 5 1 HWARNING: IMPROPER DATA ENTRY! ENTER "YES" OR "NO".) 

160 FORMAT ( 5 X ,5 OH ENTER THE RCW NUMBER OF THE ELEMENT TO BE CHANGED.) 

170 FOBMAT (5X,53HENTEE THE COLUMN NUMBER OF THE ELEMENT TO BE CHANGED 

1 ilc 
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C 



10 
20 
C — 
30 



40 



50 

60 



70 

30 

90 



C 

100 

110 

120 



130 

140 

150 

160 



SUBROUTINE READF3 (NC,NS,F3GC) 

ENTERS IHE_"C" MATRIX FEEDBACK GAI N CONTROL MATRIX 

REAL*3 FEGC ( N C , N S) ,CUM,ANSR 
INTEGER IAN3,I ,J r K,L 
DATA IY/'Y'/ f 1Z/ 1 N * / 

WRITE (5.110) 

DC 20 1= 1 , NC 
DC 10 j=i,:js 
WRITE (5,100) I,J 
CALL RDREAL (ANS3) 

FBGC (I , J) =AN S R 

CONTINUE 

CONTINUE 

CALL FSTCMS ( * CLR SCE N •) 

WRITE <5,120) 

CALL MAT PET (FBGC, NC , NS) 

WRITE (5,130) 

CALL RDCHAR (IANS) 

IF ( (IANS. NE.IY) .AND. (IANS. NE.IZ) ) GC TC 50 
GC TO 60 
WRITE (5,140) 

GC TO 40 
CONTINUE 

IF (IANS.E0.IZ) GC TC 90 
WRITE (5,150) 

CALL R DINT (IANS) 

R=IANS 

WHITE (5,160) 

CALL HDIHT (IANS) 

L= IANS 

WRITE (5,100) K , L 
CALL RDREAL ( A NS 3) 

DUM=ANSR 
DC 30 1=1, NC 
DO 70 J=1 , N S 

IF ( (I. EC-K) . AND. (J. EQ. L) ) F3 GC (I , J) =D U W 

CONTINUE 

CONTINUE 

GO TO 30 

CONTINUE 

CALL FRTCMS ('CLRSCRN «) 

RETURN 



FORMAT (5X t 14HTEE ELEMENT C ( , 1 2 f 1 H , , 12 , 2H) =) 

FORMAT I / ,oX , 5 2H ENTER THE FEEDBACK GAIN CONTROL MATRIX "C "—MATRIX 
1 .,//, 10X,44HDIMESSICN = ? CONTROLS NC X ? STATES NS .) 

FORMAT (// , 1 0 X ,45HTH E FEEDBACK GAIN CONTROL MATRIX "C"-MATRIX , // 

^ FO BHAT 1//5X, 54HDC YCtJ WISH TC CHANG £ THE VALU2 OF ANY 2AT3IX ELZH 
1ENT** // irv iQhtvctt "Vpc: » nt? 

FORMAT 
FORMAT 
FORMAT 

1 il(D 



IX, 5 1HWARNING: IMPROPER DATA ENTRY! ENTER "YES" OR "NO".) 
5 X ,5 OH E N TER THE RCW NUMBER OF TBE ELEMENT TO BE CHANGED.) 
5X.5 3 HE NT EH THE COLUMN NUMBER OF THE ELEMENT TO BE CHANGED 
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C SUBROUTINE RE A CAY *{NC,AY) 

ENTERS THE " A 11 MATRIX DIAGONAL OUT SOT COST A AT R IX . 



10 
20 
C— 
3 C 



40 

50 

60 



HEALTH AY ( NO , NC) , DUM , A 
INTEGER I AN S , I ,0 , ? , L 
DATA IY/ , Y , /,IZ/'N*/ 
WRITE (5,110) 

DC 20 1= 1 , NO 
DC 10 J= 1 , NO 
WRITS (5, ICO) I,J 
CALL RDREAL { A NS 3) 

AY (I, J)*ANSR 

CONTINUE 

CONTINUE 



Nsa 



CALL FPTCMS ( 'CLRSC5N » ) 

WRITS (5,120) 

CALL M AT PET ( A Y , N C , N C) 

WRITS (5,130) 

CALL RDCHA3 (IANS) 

IE ( (TANS. HE. II) .AND. (IANS. NE.IZ) ) GC T 
GC TO 60 
WRITS (5,140) 

GC TO 40 
CONTIS US 

IE (IANS . EO. IZ) GC TO 90 
WRITE (5.150) 

CALL 3DIST (IANS) 

K = I A N S 

WRITS (5,160) 

CALL R D I NT (IANS) 

L= IANS 

WRITE (5,100) K , L 
CALL RDREAL (ANSR) 

DUM=ANSR 
DC 30 1=1, NO 
DO 70 J = 1 , N 0 

‘ X) . AND. (J. ZQ. L) ) AY(I,J)=DUM 



C 50 



70 

80 

90 


I? ((I.: 
CONTINU1 

ccntinu; 

GO TO 3 < 
CCNTINU 
CALL F8! 
RETURN 


^00 


FORMAT 


1 10 


FORMAT 


120 


IRIX .,/ 
2 ) 

FORMAT 


130 


1. ,//) 
FORMAT 


140 


1 ENT? , // 
FORMAT 


150 


FORMAT 


160 


FORMAT 




1 .) 
END 



( 5 X , 14HTHE ELEMENT A(,I2,1I!,,I2,2H) •=) 
(//,5X ,54RENTER the outtut measurement 



COST MATRIX "A"-MAT 
X J OBSERVATIONS NO 



, 5X,53 HDIMENSION = t OBSERVATIONS NO 

(//,5X ,50HTHE OUTPUT MEASUREMENT COST MATRIX "A"-MATRIX .. 

( //5 X , 5 4 EDC YCU WISH TO CHANGE THE VALUE 0? ANY MATRIX ELZM 
, 1 CX . 1 9HT Y PE "YES' 1 OR "NO",) 

( IX, 5 1 HWA5NING: IMPROPER DATA ENTRY! ENTER "YES" OR "NO".) 
5X,50HENTER THE RCW NUMBER OF THE ELEMENT TO BE CHANGED. [ 
5X,53HENTER THE COLUMN NUMBER OF THE ELEMENT TO 3E CHANGED 
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c 

c== 



10 
c — 
2 0 



30 

40 

50 



60 

70 



SUBROUTINE 2EAD3 (NC,3) 

ENTERS THE " 3 " MAT HI X CO NT R C I^COST WEI GHT ING_ MAT RIX 

RFAL*3 £ (NC , N C) , CUM , AN r i R 
INTEGER IABS r I # J # ? f L 
DATA I Y/ 1 Y 1 / , 12/* N f / 

WHITE (5,90) 

DC 10 1=1, NC 
DC 10 J= 1 ,NC 
W3IIZ (5 , cO) I ,J 
CALL 3DHEAL (ANS3) 

B (I,J) =A NSR 



C 

80 

90 

100 

110 

120 

130 

140 



CALL FRTCMS ( * CLE SC E N ') 

WHITE (5 f 110) 

CALL LATEST ( E ,NC,NC) 

SHITE (5,110) 

CALL RDCHAR (IANS) 

I? ((IANS.NE.il) .AND. (IANS. NE* 12)) GC TC 40 
GC TO 50 
SHITE (5 ,12 0) 

GC TO 30 
CONTIN US 

IF (IANS.EQ.IZ) GC TC 70 
SHITE (5,130) 

CALL RDINT (IANS) 

K = I A N S 

SHITE (5,140) 

CALL 3DINT (IANS) 

L = IANS 

WHITE (5, SO) X,L 
CALL* SDHFAL ( A NS R) 

DUM=ANSa 
DC 60 1=1, NC 
DO 60 J = 1 , N C 

IF ((I.EC-K) -ANC. (J. EQ. L) ) 3 (I, J) =00 2 
CONTIN UE 
GC TO 20 
CONTINUE 

CALL FRTCMS ( 9 CL ESCB N ») 

RETURN 



FORMAT (5X,14HTHZ ELEMENT 3 ( , 12 , 1 H , . 12 , 2H) =) 

F 0 fi £ AT (/ , b X , 52EENTF? THE CONTROL COST WEIGHTING MATRIX "3"-2ATRI 
IX ,/, 1 QX,45HDIMENSICN = * CONTrCLS NC X ♦ CONTROLS NC ) 

FORMAT (// .1 0X.37HTHE CONTROL CCST M AT E IX 3...,//) 

FOEMAT (//5X , 54 EDC YOU WISH TC CHANGE THE VALUE OF ANY MATRIX ELEM 



12NT?,//, 1CX, 1 9HTYFS "YES " OR '’NC’*.) 

( IX, 5 1H WARNING: IMPROPER DATA ENTRY! 



fcem&i 

FORMAT 

FCEMAT 

1 ^ND 



ENTER "YES” OR "NO" . ) 



(5X,50HENTFR IHE 2 CW NUMBER OF THE ELEMENT TO 3E CHANGED.) 
(5X, 5 2 B ENT FR THE COLUMN NUM3ER OF THE ELEMENT TO 3E CHANGED 
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SUBROUTINE READG2 (NS,N G , IG Art , GArt) 

ENTERS THE "GAM" MATRIX IROCESS NOISE DISTR I3UTI0N 3 ATS IX . 



10 
20 
C — 
3 0 
40 



50 



60 

70 



80 

90 

100 



C— ■ 
110 
120 



130 

140 



150 

160 

170 



RE AL* 3 GArt (NS, 21 G) ,CCrt,ANSR 



INTEGER IANS,I,J,K,L 
DATA n/ r <'/ { it/ 1 N * / 
IF JIG AH . SQ. 1) GC TO 
WRITE (5,120) 

DC 20 1=1, MS 
DO 10 J= 1 f NG 
WRITE (5,110) I,J 
CALL BOREAL ( A NS R ) 

G AH (I, J) - ANS 3 
CONTIU US 
CONTINUE 



) Art 



40 



CALL FRTCMS (*CL3SCHN •) 

CONTIN 02 
WHITE (5,130) 

CALL rtATFRT (G Art , NS , NG) 

WHITE (5,140) 

CALL BDCHAR (IANS) 

IF ( (IANS. NE. IY) .AND. (IANS. ME. IZ) ) GO TO 60 
GO TO 70 
WHITE (5,150) 

GO ro 50 
CONTINUE 

IF (IANS. 20. IZ) GC TC 100 
WHITE (5,160) 

CALL EDINT (IANS) 

K=I ANS 

WRITE (5,170) 

CALL R D I NT (IANS) 

L = I A N S 

WHITE (5,110) K , L 
CALL RDREAL (AKSR) 

DU H= ANS B 
DO 90 1=1, NS 
DC 80 J = 1 , N G 

I? Ml • E Q • K) -AND. (J. EQ. L) ) GArt (I, J) = DUrt 

CONTINUE 

CON TIN US 

GC TO 30 

CONTINUE 

CALL F ST CMS ('CLRSCHN ' ) 

RETURN 



FORMAT (5X t 16HTHE ELEMENT G Art ( ,12 , 1H , , 12, 2 H) =) 

F03MAT /,5X, 36EESTFR THE PROCESS NOISE D ISTRI 3UTIO N ,/ , 5X , 2 4 H 25 AT RI 

IX ”G A M M A M A T B I X . ,/, 2X , 56 HDIMEUSIO N = * STATES NS X # PROCESS 

2NOISE SO UHCE S NG ) 

FORMAT (//,10X,37HTHS PROCESS SCISE DISTRIBUTION MATRIX,/, 10X,19H 
1 11 GArt - MATRIX //) 

FORMAT (//5 X, 5 4HQC ’ Y CD WISH TC CHANGE THE VALUE OF ANY HATH IX SLEH 
1ENT?,//, 1CX, 19ETYEE "YES" OH "NO".) 

FOHHAT (IX, 51 RViAHNING: I HF80 ?E3 DATA ENTHY! ENTER "YES" OH "NO".) 
F OHM AT (5X,50E£NTER THE HOW NUHBER OF THE ELEMENT TO BE CHANGED.) 

FOHHAT (5X,53HENTEH THE CCLUHN NIIHBE3 CF THE ELEMENT TO 3E CHANGED 

1 .) 



END 
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SUBROUTINE REATQ |NG,Q) 

i5 E _" Q l!_ XATEI ?L ^2 TEE *eightisg_:i atri x = 

RE AL*8 C(NG,NC*) ,EUM,ANSR 
INTEGER IANS,I,J,R,L 
DATA I Y/ • Y 1 / , IZ/ f N V 
WRITS (5.110) 

DC 20 1=1, NG 
DC 10 J=1,NG 
WRITS (5,100) I,J 
CALL RDREAL ( A NS3) 

Q <I,J) = ANSR 
CONTINUE 
CONTINUE 



10 
20 
C — 
3 G 



40 

50 

60 



70 

30 

90 



CALL FRTCMS ( * CLRSCS N * ) 

WHITE (5 , 120) 

CALL MATFRT ( Q ,NG ,NG) 

WRITE (5,130) 

CALL R DC HAS (IANS) 

I? jJ (IANS. ME. IY) .AND. (IANS. HE.IZ) ) GC TC 50 
G C i 0 60 
WRITE (5,140) 

GC TO 40 
CONTIN UE 

IF (IANS . SO. IZ) GC TO 90 
WRITE (5,150) 

CALL R DINT (IANS) 

K = I A N S 

WHITE (5,160) 

CALL RDINT (IANS) 

L- IANS 

WRITS (5,100) K,L 
CALL iDREAL ( A NS R) 

DUM=ANSR 
DO 80 1=1, MG 
DO 70 J=1 , NG 

IF ( (I. 2C-X) • AND. (J. cQ. L) ) Q (I ,J) =DU M 

CONTINUE 

CONTINUE 

GO TO 3 0 

CONTINUE 

CALL FPTCMS ( 1 CIS SC5 N ') 

BETURN 



C — 

10 0 FORMAT (5X, 1 4 HT F E ELEMENT Q ( , 1 2 , 1 H , , 12 , 2H) =) 

110 FORMAT ( / / , 5 X , U 4 H E N T EE THE PROCESS NOISE PSD WEIGHTING MATRIX,/, 5X 

1 , 1 2H " 0 " M ATR I X . ,// , 5X, 4 2 HD I ME bSION = * PROCESS NOISE SOURCES NG 
2 X,/,1 7X ,27H 4FFCCISS NOISE SOURCES NG ) 

120 FORMAT ( / / , 5 X , 4 2 HTHS PROCESS SCISE WEIGHTING MATRIX Q..,//) 

130 FORMAT (//5X,$4HCC YOU WISH TO CHANGE THE VALUE OF ANY MATRIX ELEM 
1SNT?, //, 10X, 1 9HTYFE "YES" OR "NC".) 

140 FORMAT ( 1X,5 1HWAHNING: IMPROPER DATA ENTRY! ENTER "YES" OR "NO"-) 

150 FORMAT 5X,50HENTER THE RCW NUMBER OF THE ELEMENT TO 3E CHANGED.) 

160 FORMAT (5X,53HENTER THE COLUMN NUMBER OF THE ELEMENT TO BE CHANGED 

1 E^D 
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subroutine reach {ho ,rc) 

ENTERS THE " R " MATRIX_ HEaElTT^ NOXSE__ DISTRIBJTION MATRIX .. = 

BEALES EC ( NO , NO) ,DUM,ANSR 



10 
C — 
2 C 



30 

40 

50 



INTEGER IANS, 



H f f 



DATA n/»YV, 

WRITS (S f SC) 

DC 10 1=1, NO 
DC 10 J = 1 1 N G 
WRITE (5 ,80) I.o 
CALL RDREAL ( A NSR) 
RC (I, J) =ANSH 



K,L 



/‘SV 



60 

70 



CALL ? RTCMS ( • CLHSCS N ') 
WRITE (5,100) 

CALL MATPRT (RC,NC,NC) 
WRITE (5,110) 

CALL RDCHAR (IANS) 



C — 

30 

90 



100 

110 

120 

130 

140 



I? J (I ANS.N2 . IX) • AND . (IANS. N£. IZ) ) GC TC 40 
GC TO 50 
WRITE (5,120) 

GC TO 30 
CONTIN 02 

IF (I A NS • FQ . IZ ) GC TC 70 
WRITE (5,130) 

CALL R DINT (IANS) 

K=IANS 

WRITE (5,140) 

CALL RDINl (IANS) 

L= IANS 

WRITE (5,80) K,L 
CALL 3 DR £ AL ( A NSR) 

DUM=ANSR 
DC 60 1=1, NO 
DO 60 J = 1 , N O 

IF ( (I. EQ.X) • AND. (J. SQ. L) ) RC {I, J) =DUM 

GO TO 20 

CONTINUE 

CALL ? RTCMS ( 1 CIR SC 5 N «) 

RETURN 



FORMAT (5X t 14RTEc ELEMENT R(,I2,1H,,I2-2H)=) 

FORMAT f/,Sx, 6CHHNTER THE MEASUREMENT NCISE DISTRIBUTION MATRIX " 
1 R M M ATR I X • ,// ,5X,53HDIMENSION = * OBSERVATIONS NO X t OBSERVATIO 
2NS NO ) 

FORMAT ( // , 1 5 X ,5 CHTHE MEASUREMENT NOISE DISTRIBUTION MATRIX R. 

FORMAT (//5X.5UHDC YCU WISH TC CHANGE THE VALUE OF ANY MATRIX EL EH 

IENT?,//, 1CX, 19H1YEE "YES " OR "NO".) 

... ... 
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SUBROUTINE REALFE (NS,N0,?3G2) 

I NTEHACT I VELY ENTERS THE "K " FEEDBACK GAIN ESTIMATOR MATRIX 



3 C CALL PRICES (‘CLHSCIN *) 

WRITE (5,120) 

CALL MATEBT ( F EG E , NS , NO) 

40 WRITS (5.130) 

CALL R DC fc AR (IANS) 

IF ((IAMS.arS.il) .AND. (IANS. NE.IZ) ) GC TC 50 
GC TO 60 

50 WRITE (5,140) 

GC TO 40 

60 CONTINUE 

IF (IANS . EQ. 12) GC TC 90 
WRITS (5.150) 

CALL SEIST (IANS) 

K = I A N S 

WRITS (5,160) 

CALL RDINT (IANS) 

L = IANS 

WRITE (5,100) *,L 

CALL RDREAL (ANSE) 

D UM = A N SR 
DC 30 1=1, NS 
DO 70 J - 1 , N 0 

IF ( (I. EC-K) . ANE. (J. EO.L) ) F B GE ( I , J) =D 0 tt 
70 CONTINUE 

30 CONTINUE 

GO TO 30 
90 CONTINUE 

CALL ? RTCMS (’CLHSCEN ’) 

RETURN 



100 FORMAT (5X.14HTEE ELEMENT K ( , 1 2 , 1 H , , 12 , 2H) =) 

110 POSMAT ( / t 5X , 54HENTER THE FElLEACK GAIN ESTIMATOR MATRIX "K" -MATS 
1IX ..//, 1CX. 48HDIKENSI0N = 4 STATES NS X 4 OBSERVATIONS NO .) 
120 FORMAT ( // , 1 5 X , 4 7HTH E FEEDBACK GAIN ESTIMATOR MATRIX "K"- MATRIX , 
^//) 

130 FORMAT (//, 5 X , 5 4 H CO TOU WISH TC CHANGE THE VALUE OF ANY MATRIX EIE 
1 MENT? , // . 10X. 1 9HTKPE '•YES** CR "NO". ) 

140 FORMAT ( IX, 5 1 HWARNING: IMPROPER DATA ENTRY! ENTER "YES" OR "NC".) 

150 FORMAT ( 5 X ,5 0 HENT EH THE R CW NUMBER OF THE ELEMENT TO BE CHANGED.) 

'T'n p r rHifJCFn 




DC 10 J— 1 , NO 
WRITE (5.100) I,J 
CALL RDR ill ( A NS R ) 



NC) , CUM, A NS R 



iu 

20 U « i. i. .1 U C. 

C-* 
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c 



10 
c- 
2 C 

30 

40 

50 



60 

70 



C — ■ 
80 
90 
100 

110 



120 

130 

140 



SUBROUTINE REAEW <NG,WR) 

I NTERACT IVEL Y ENTERS "HO" STEADY DISTURBANCE VECTOR 



SEAL*3 WF(NG) ,CUM,AKSR 
INTEGER IANS, I. K 
DATA lY/’YV, 12/ * N * / 
WRITE (5.100) 

DC 10 1-1. MG 
WRITE (5,30) I 
CALL 3DREAL ( A N S 2 ) 

W H (I) = AN SB 
CONTIN UE 



CALL F R T C M S ( 1 CI3SC5 N • ) 

WHITE (5,110) 

WRITE 5,90) (WR (I) ,1=1 ,NG) 

WRITE (5.12 0) 

CALL RCChAR (IANS) 

IF ^lANS.NE.IY) .AND. (IANS. NE.IZ) ) GC TC 40 

WRITE (5,130) 

GC TO 30 
CONTIN UE 

IF (IANS. EQ. 12) GC TC 70 
WRITS (5,140) 

CALL RDINT (IANS) 

K = 1 AN S 

SHITE (5,80) K 
CALL RDREAL (AKSP) 

DUM= ANSR 

DO 60 1*1, MG 

IF (I. EQ.S) W9(I)=CUM 

CONTIN UE 

GC TO 20 

CONTINUE 

CALL FETCMS ('CLRSC3N •) 

RETURN 



FORMAT (5X.15HTHE ELEMENT W0(,I2,2H)=) 

FORMAT { E 12. 5 ) 

FORMAT (/,5X,57HENTER THE STEADY DISTUREANCE VECTOR MATRIX "W0"-M 
1ATRIX . ,//,l6x,44HDIMENSICN = * PROCESS NOISE SOURCES NG X 1) 
FORMA! ( //,1 5 X,53HTHE STEADY DISTURBANCE VECTOR MATRIX "V0"-«AT3I 

IX ...,// 

FORMAT ( //5 X , 5 4 H DC YCU WISH TC CHANGE THE V ALU E OF ANY MATRIX EL EM 
1 ENT?,//, 10X, 19RTYPE "YES " OR '•NO".) 

FORMAT ( IX, ol HWARNING: IMPROPER DATA ENTRY! ENTER "YES" OR "NO".) 
FORMAT ( 5X,5 CHENIER I HE ROW NUK3E3 OF I HE ELEMENT TO BE CHANGED.) 
END 
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SUBROUTINE RDBEAL — INTER ACTIVELY BEADS A BEAL NUH5ER REPLY 

INTO A FORTRAN PROGRAM. IF THE USER IN AD V E3TEN LY ENTERS A NULL = 

SIRING THE 5/ H ISSUES A WARNING AND ALLOWS A RECOVERY, 

SUBROUTINE RDREAL ( A AS 3) 

EEAL*8 ANSR 
INTEGER COUNT 

CO UNT= 0 

10 CONTINUE 

CQONT-COUNT+1 
I? (COUNT.LT. 3) GC TC 20 
WRITS (5,60) 

GO TO 40 
20 CONTINUE 

READ (5,*,END=3Q,E5R=3Q) ANSR 
RETURN 

30 REWIND 5 

WRITE (5,50) 

GC TO 10 
40 CONTINUE 

STOP 

50 FORMAT ( IX, 64 HWARNING: NULL STRINGS ARE NOT ALLOWED, ENTER A NU HE 

13ICAL VALUE.) 

60 FORMAT <////, 2X, 4 2HP HOG RAM KILLED - TWO NOLL STRINGS ENTERED!,/) 

END 
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SUBROUTINE 3 DINT — INTERACTIVELY PEADS AN INTEGER REPLY 

INTO A FORTRAN PROGRAM- IF THE USER I N ADV ERTSNLY ENTERS AN IMPROPER= 

DATA CHARACTER THE S/R ISSUES A EARNING AND ALLOWS A RECOVERY. 



SUBROUTINE R D INT (IANS) 
INTEGER COUNT , IANS 




COUNT=0 






10 


CONTINUE 

CCCNT^COUNT* 1 

IF (COUNT.LT. 3) GO TC 20 








WRITS (5 ,60) 
GO TO 50 






20 


CONTIN UE 

READ (5, *,END=40, £RR=40) 
I? (IANS) 40,40,30 


IANS 




30 


CONTINUE 

RETURN 






40 


REWIND 5 
WRITE (5,70) 
GO TO 10 






50 


CONTINUE 

STOP 






60 


PCHaAT (//,5X,495EECGaAa 


TERMINATION' 


- TWO IMPROPER DATA ENTRIES!! 


70 


1 1 

FORMAT ( IX , 5 6H W A3NING: I 
1GE3.) 


M PRO PER DATA 


ENTRY! ENTER A POSITIVE INTE 




END 
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SUEHCUTINS BLZ'd A S — INTER ACTI7 ELY READS A CHARACTER STRING REPLY - 

('YES' OR ' NC») I STC A FORTRAN F 5CG3AK. IF THE USER IN ADV ERTENL Y 
EST2BS A HULL STBI2IG_TEE_S/H ISSUES A EARNING AND ALLOWS \ R SCO VER Y= 

SUBROUTINE RDCHAE (IANS) 

INTEGER CCU NT , IANS 
DATA X Y/ ' I 1 / , I Z / ' S'/ 



10 



20 



30 



40 

C 

50 

60 

70 



CCCNT=0 

CONTINUE 

COUNT=CO CNT+ 1 

IF (COUNT. LT. 3) GC TC 20 

WRITE (S r 60) 

GO TO 40 
CONTIN UE 
REWIND 5 

READ (5,7C,ZND=30,2R R-3 0 ) IANS 

RETURN 

REWIND 5 

WRITE (5,50) 

GC TO TO 
CONTINUE 
STOP 



FORMAT (1X,6QH WARNING; NULL STRINGS ARE NOT ALLOWED, ENTER "YES" 

1 C 5 "SO 11 • ) 

FORMAT ’(///, 5X,U7HPRCGRAM TERMINATION - TWO NULL STRINGS ENTERED!) 

T (A if 



FORMAT 
END 



noon 



SUBROUTINE SATPRT -- DISPLAYS A TSC-DI M E NS ION AL ARRAY (16 COLS. MAX) = 
IN 7 ABIABL2 SCREEN FORMAT FOR USEE EASE IN ROM IDENTIFICATION. = 



SUBROUTINE MATEET ( E BTT , N SOM , SCCL) 
IMPLICIT RE AL * 8 (A-B,0-2) 

DIMENSION ?P. IT (iJFOW ,NCOL) 



IP ( 


NCOL.EQ.0) 


1 


NCCL- 1 












IF 1 


[NCOL.EQ. 1] 


1 


WHITE ( 


5, 10) 


( (PRTT (I, J) 


, J = 1 , MCO L) 


,1=1 


,NROW) 


IF { 


[NCOL.EQ. 2] 


1 


SHITE ( 


5,20) 


( ( ? RT T (I, J) 


,J=1, ^COL) 


,1=1 


, N R 0 W ) 


IP { 


NCOL . EQ . 3 ] 


1 


WRITS 


5 ,30} 


( JPSPT (I, J) 


,j- i , :icoi) 


,1=1 


,K30W{ 


I? 1 


NCOL. EQ. 4 1 


1 


WHITE ( 


'5, U0( 


i (P?T? (I, J) 


,J=1 , NCOL) 


,1=1 


, M R 0 a ) 


IF i 


ncol.eo.5; 




WRITE { 


[5 ,50) 


< (PRTT (I , J) 


, j= i, :icol) 


, 1=1 


, N R 0 W ) 


IF < 


[NCOL. HO. 6 ] 




WRITS ( 


5, 60) 


( ?HTT (I,J) 


,j= i,:icod 


,1= 1 


, N R 0 W ) 


IP \ 


[NCOL.EQ.7 




WHITE ( 


5,70) 


( (PRTT I, J) 


,J*1 ,MCOL) 


,1=1 


,NROW) 


IP { 


[NCOL . EQ . 3 




WHITE 1 


[5 ,30) 


( (PRTT (I , J) 


,J=1, NCOL) 


,1=1 


, N H 0 W ) 


IF I 


[NCOL. EQ. 9[ 


1 


WHITE ( 


5,90) 


( (PRTT (I, J) 


,J= 1 , WOOL) 


,1=1 


, NRO W) 


IF 1 


[NCOL.EO. 1 ■ 


01 


) WHITE ( 


5 ,100) 


( (PRTT (I , J) 


, J- 1 , NCO Lj 


,1=1 


, :jrow) 


IP 1 


[NCOL. EQ. 1 1 [ 


> WHITE i 


5, 110) 


( (PRTT (I, J) 


,J=1, NCOL 


,1=1 


, SHOW) 


IF 1 


[NCOL. EQ. 12] 


} WHITE ( 


[5, 120) 


( (PRTT (I, J) 


, J- 1 , NCOL) 


,1=1 


, il R 0 W ) 


IF | 


(NCOL.EQ. 1 3 


) WHITE i 


(5,130) 


( (PRTT {I , J) 


, Jsl , NCOL) 


,1=1 


,'JROW) 


IF I 


NCOL. EO. 14 


) W B IT E \ 


[ 5 , 14 0) 


((PRTT (I,J) 


, J=1 , RCOL) 


,1=1 


,NROW) 


IF i 


NCOL.EQ. 1 5 


) WHITE | 


[5 , 15 0) 


( (PRTT (I, J) 


,J= 1, HCOL) 


,1=1 


,NROW) 


IF i 


(NCOL.EQ. 16 


WHITE . 


(5,160) 


( (PRTT (I , J) 


, J- 1, ^COL) 


,1=1 


,NROW) 


BETU3M 

















10 


FORMAT 


(512- 5) 




20 


FORMAT 


f 2F12 . 5) 




30 


FCSM AT 


(3512- 5 




40 


FORMAT 


( 4 F 1 2 . 5 




50 


FORM AT 


(5512. 5) 




60 


FORMAT 


(6F12.5) 




70 


FCfiflAT 


(6r12.5,/,F12-5, //) 




30 


FORMAT 


(6?12.5,/,2F12.§, //} 




90 


F OHM AT 


6F12. 5,/, 2512.5,// 




100 


FORMAT 


6F12.5,/,4F12.5 , // 




110 


F CBM AT 


(6512. 5,/,5F12.5, //j 




120 


FORMAT 


(6F12 - 5,/,6F1 2.5 ,//' 




130 


FORMAT 


(6F12. 5,/,6F12.5,/,l 


F 1 2 • 5 .//) 


140 


FORMAT 


6F12.5,/,6F12.5,/,: 


2F1 2 . 5 ,//) 


150 


FCSM AT 


( 6 F 1 2 . 5,/,6F'12.5,/, 


3F12.5,// 


160 


FORMAT 


<6F12.5,/,6F12.S,/,4F12.5,//) 



ESC 
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SUBROUTINE M ATPRT — DISPLAYS A TWC-DI H E NS ION AL ARRAY (16 COLS, MAX )- 
IN VARIABLE SCREEN FORMAT FOR USER EASE IN SOW IDENTIFICATION. 



SUBROUTINE MATP5T (FSTT . NROW , NCCL) 
IMPLICIT SEAL * 8 (A-H,0-Z) 

DIMENSION PRTT <NFOW,HCOL) 



C 

10 

20 

30 

40 

50 

60 

70 

30 

90 

100 

110 

120 

130 

140 

150 

160 



I? 

IF 

IF 

IF 

IF 



NCOL « 
NCOL . 
NCOL . 
NCOL. 
NCOL. 
I? (NCOL. 
IF (NCOL. 
IF (NCOL. 
IF (NCOL. 
IF (NCOL. 
IF (NCOL. 
IF (NCOL. 
““ NCOL. 



I? 

IF 

IF 

IF 

IF 



S NCOL. 
NCOL. 
NCOL. 
NCOL. 



RETURN 




WRITE ( 


5, 10) 


( (PRTT (I,J) 


WRITE ( 


5 ,20) 


( ( PSTT (I, J) 


WRITE 


;s,3o 


( (PRTT (I, J 


WRITE 


5 J 40 ( 


( (PSTT (I, J 


WRITE j 


;s,so) 


((PRTT (I, J 


WRITE ( 


5,60 


( (PRTT (I , J) 


WRITE ( 


5 , 70 ) 


( (PRTT (I, J) 


WRITE | 


[5,30) 


((PRTT (I,J) 


WRITE ( 


5,90) 


( (PRTT (I, J) 


WRITE ( 


5 ,100) 


((PRTT (I, J) 


WRITE l 


(5, 110) 


( (PRTT (I, J) 


WRITE ( 


; 5 , 12 0) 


((PRTT (I,J) 


WRITE i 


(5,130) 


( (P3TT (I, J) 


WRITE | 


[5, 14 0) 


((PRTT (I,J) 


WRITE < 


5 , 150) 


( (PRTT (I, J) 


WRITE ' 


(5 ,160) 


((PRTT ( i # j 



'1, NCOL) ,1 = 
■1, NCOL) , 1 = 

• 1 , NCOL) ,1= 
= 1 ,NCOL) ,1 = 
1, NCOL) , 1= 
: 1 , NCO L) ,1 = 
*1,NCOL) ,1= 
: 1 , NCOL) , I- 
1 1 , NCO L) ,1 = 
: 1 , NCOL) ,1= 
: 1 , NCOL) ,1= 
■1,NCO L) ,1= 
-1 , NCOL) , 1= 
= 1 , NCOL) ,1= 
: 1 , NCOL) ,1= 
= 1 , NCOL) ,1= 



1 ,NROW) 
1 , SHOW) 
1,NHOWj 
1 , NROW) 
1,NKOW) 
1 , N 3 O W ) 
1 ,NROW) 
1 , N ROW) 
1 f NROW) 
1, NROW) 
1,NROW) 
1 , N ROW) 
1,NROW) 
1 ,NROW) 
1 ,NROW) 
1 , N R O W ) 



FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

ESC 



(f 12. 
(2F12 
( 3f 1 2 
( 4 F 1 2 
( 5 F 1 2 
( 6F1 2 
6F12 
6F12 
6F12 
( 6F1 2 
( 6 F 1 2 
(6F12 
( 6 F 1 2 
(6F12 
( 6 FI 2 
( 6F1 2 



% 

• 5 
.5 
.5 

• | • / 
• 1 ’/ 

: V y 

. S,/ 



F12.5 t //) 
,2?12.5 ,// 

3F 12.5,//: 

,4?1 2.5 ,// 

, 5F 1 2. 5 , // 
r 6F 1 2,5 ,//) 

, 6F 1 2. 5 , / ,r 1 2 . 



. . . <..^,//) 
6F12.5,/,2F12.5 , //) 
6F12.5,/,3F12.5,//S 
6F12.5,/,4F1 2.5,//) 
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